Read in the data spreadsheets:
sheet <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/LiveFishInteraction_Feb9.csv", sep = ",", strip.white = TRUE)
scored <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/ScoredFishInteraction_Feb19.csv", sep = ",", strip.white = TRUE)
# Combining the two is easier when dealing with the relevant columns alone. These are:
# Species.A, Species.B, Group.A, Group.B, and Winner.
ABwin_live <- sheet[,c(4,6,8,9,10)]
ABwin_scored <- scored[,c(10,12,13,15,16)]
both <- rbind(ABwin_live, ABwin_scored)
# Create new versions of the three data frames that include one-on-one interactions only:
sheet_oneonone <- sheet[sheet$Group.A == 1 & sheet$Group.B == 1,]
scored_oneonone <- scored[scored$Group.A == 1 & scored$Group.B == 1,]
both_oneonone <- both[both$Group.A == 1 & both$Group.B == 1,]
Step 1. Unidentified species are marked with an asterisk. Exclude them from the matrix based on this criterion:
drop_ufos <- function(df) {
# Dummy vector for storing the numbers of the rows with unidentifiable species:
tobedropped <- vector()
for (i in 1:nrow(df)) {
if (grepl("*", df$Species.A[i], fixed = TRUE) == TRUE |
grepl("*", df$Species.B[i], fixed = TRUE) == TRUE |
grepl("*", df$Winner[i], fixed = TRUE) == TRUE) {
tobedropped <- c(tobedropped, i)
}
}
return(df[-c(tobedropped),])
}
sheet2 <- drop_ufos(sheet)
scored2 <- drop_ufos(scored)
both2 <- drop_ufos(both)
sheet_ooo2 <- drop_ufos(sheet_oneonone)
scored_ooo2 <- drop_ufos(scored_oneonone)
both_ooo2 <- drop_ufos(both_oneonone)
Step 2. Remove intraspecific interactions (a few of these were scored for live observations in the beginning, before it was decided to selectively focus on interspecific interactions only):
drop_intra <- function(df) {
tobedropped <- vector()
for (i in 1:nrow(df)) {
if (identical(as.character(df$Species.A[i]), as.character(df$Species.B[i])) == TRUE) { tobedropped <- c(tobedropped, i)
}
}
if (length(tobedropped) > 0) {
return(df[-c(tobedropped),])
} else {
return(df)
}
}
sheet3 <- drop_intra(sheet2)
scored3 <- drop_intra(scored2)
both3 <- drop_intra(both2)
sheet_ooo3 <- drop_intra(sheet_ooo2)
scored_ooo3 <- drop_intra(scored_ooo2)
both_ooo3 <- drop_intra(both_ooo2)
Step 3. Drop “Striated Surgeonfish” from the matrices, as this is likely a composite category comprising both striated and brown surgeonfish. These species may have different dominance ranks, but are almost indistinguishable in appearance.
no_surgeon <- function(df) {
return(df[df$Species.A != "Striated Surgeonfish" & df$Species.B != "Striated Surgeonfish",])
}
sheet4 <- no_surgeon(sheet3)
scored4 <- no_surgeon(scored3)
both4 <- no_surgeon(both3)
sheet_ooo4 <- no_surgeon(sheet_ooo3)
scored_ooo4 <- no_surgeon(scored_ooo3)
both_ooo4 <- no_surgeon(both_ooo3)
Extract the levels of Species.A and Species.B and use them as the rows and columns of the dominance matrix:
# The union of the species listed in the "Species A" and "Species B" columns will be used
# as the rows and columns of a new matrix filled with zeros:
emptymatrix <- function(df) {
species <- union(levels(factor(df$Species.A)), levels(factor(df$Species.B)))
rownum <- length(species)
newmatrix <- matrix(0, nrow = rownum, ncol = rownum)
row.names(newmatrix) <- species
colnames(newmatrix) <- species
return(newmatrix)
}
dommatrix <- emptymatrix(sheet4)
scoredmatrix <- emptymatrix(scored4)
bothmatrix <- emptymatrix(both4)
domooo <- emptymatrix(sheet_ooo4)
scoredooo <- emptymatrix(scored_ooo4)
bothooo <- emptymatrix(both_ooo4)
# If species A won over species B (winner = A), add 1 to the cell in column A and row B.
# If species B won over species A (winner = B), add 1 to the cell in column B and row A.
# An increment function will be used to increase the value of each cell accordingly:
inc <- function(x) # credit: Grega Kešpret, Stack Overflow
{
eval.parent(substitute(x <- x + 1))
}
# Now apply it to populate the matrix:
populate <- function(sheet, df) {
for (i in 1:nrow(sheet)) {
spec_a <- as.character(sheet$Species.A[i])
spec_b <- as.character(sheet$Species.B[i])
winner <- as.character(sheet$Winner[i])
if (identical(spec_a, winner) == TRUE) {
inc(df[as.character(spec_b), as.character(spec_a)])
}
if (identical(spec_b, winner) == TRUE) { # This should be equivalent to "else"
inc(df[as.character(spec_a), as.character(spec_b)])
}
}
diag(df) <- NA
return(df)
}
dommatrix <- populate(sheet4, dommatrix)
scoredmatrix <- populate(scored4, scoredmatrix)
bothmatrix <- populate(both4, bothmatrix)
domooo <- populate(sheet_ooo4, domooo)
scoredooo <- populate(scored_ooo4, scoredooo)
bothooo <- populate(both_ooo4, bothooo)
Convert from matrices to data frames for easier manipulation:
convert_to_df <- function(df) {
frame <- as.data.frame(df, stringsAsFactors = FALSE)
for (i in 1:ncol(df)) {
frame[,i] <- as.numeric(frame[,i])
}
return(frame)
}
domall <- convert_to_df(dommatrix)
scoredall <- convert_to_df(scoredmatrix)
bothall <- convert_to_df(bothmatrix)
liveoooall <- convert_to_df(domooo)
scoredoooall <- convert_to_df(scoredooo)
bothoooall <- convert_to_df(bothooo)
Exclude cornetfish, the only predator in the dataset. (Note: the fact that this is done here, at the matrix formatting stage rather than at the data clean-up stage, simply reflects the point at which the decision was taken during the actual analysis. Excluding a taxon from the matrix using the method below is equivalent to excluding it from the data sheet.)
# Exclude a taxon from a matrix:
# Original unvectorized function used in the pre-revisions part of the file is commented
# out below:
#
# exclude <- function(df, taxon) {
# colrow <- c(which(row.names(df) == taxon), which(colnames(df) == taxon))
# excluded <- df[-colrow[1], -colrow[2]]
# return(excluded)
# }
exclude <- function(df, taxon) {
retained <- df[!rownames(df) %in% taxon, !colnames(df) %in% taxon]
return(retained)
}
# Note: since there are no interactions involving cornetfish in the video-only data,
# we only need to exclude it from the live-only and combined datasets.
nocornet_domall <- exclude(domall, "Cornetfish")
nocornet_bothall <- exclude(bothall, "Cornetfish")
nocornet_live <- exclude(liveoooall, "Cornetfish")
nocornet_both <- exclude(bothoooall, "Cornetfish")
compl <- function(df) {
counter <- 0
for (i in 1:nrow(df)) {
for (j in 1:ncol(df)) {
if (!is.na(df[i,j]) & df[i,j] == 0 & df[j,i] == 0) {
counter <- counter + 1
}
}
}
return(100*(1 - (counter/(nrow(df)*(nrow(df) - 1)))))
}
paste("The completeness of the live-only all-species matrix is ", compl(nocornet_domall), "%", sep = ""); paste("The completeness of the video-only all-species matrix is ", compl(scoredall), "%", sep = ""); paste("The completeness of the combined all-species matrix is ", compl(nocornet_bothall), "%", sep = ""); paste("The completeness of the live-only one-on-one all-species matrix is ", compl(nocornet_live), "%", sep = ""); paste("The completeness of the video-only one-on-one all-species matrix is ", compl(scoredoooall), "%", sep = ""); paste("The completeness of the combined one-on-one all-species matrix is ", compl(nocornet_both), "%", sep = "")
## [1] "The completeness of the live-only all-species matrix is 38.961038961039%"
## [1] "The completeness of the video-only all-species matrix is 27.2058823529412%"
## [1] "The completeness of the combined all-species matrix is 37.9446640316206%"
## [1] "The completeness of the live-only one-on-one all-species matrix is 33.7662337662338%"
## [1] "The completeness of the video-only one-on-one all-species matrix is 25%"
## [1] "The completeness of the combined one-on-one all-species matrix is 33.596837944664%"
library(gdata)
library(compete)
# To determine linearity, we need to make our dominance matrix symmetric. Below, we do
# this by subtracting the lower triangle from the upper triangle, and flipping the results
# around the diagonal:
symmatrix <- function(df) {
# Going through the upper triangle column-wise corresponds to going through the lower
# triangle row-wise:
triangle <- as.numeric(upperTriangle(df)) - as.numeric(lowerTriangle(df), byrow = TRUE)
# Replace the upper triangle in the default column-wise order:
symmat <- matrix(NA, nrow = nrow(df), ncol = ncol(df))
upperTriangle(symmat) <- triangle
# Replace the lower triangle in the row-wise order:
lowerTriangle(symmat, byrow = TRUE) <- triangle
return(symmat)
}
# Note that the p-values are simulation-based, and therefore stochastic: they will differ
# slightly between runs
devries(nocornet_domall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.3196588
## p-value= 0.0015
devries(scoredall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2982113
## p-value= 0.0528
devries(nocornet_bothall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2962757
## p-value= 0.002
devries(nocornet_live, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2779799
## p-value= 0.0097
devries(scoredoooall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2972647
## p-value= 0.0559
devries(nocornet_both, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2678731
## p-value= 0.0066
Uncomment the code below (and change the paths as needed) to print the matrices:
# write.table(nocornet_domall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/live_only_matrix.txt", quote = FALSE, sep = "\t")
# write.table(scoredall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/video_only_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_bothall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/both_datasets_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_live, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/live_only_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
# write.table(scoredoooall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/video_only_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_both, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/both_datasets_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
The code below was used to compute the Clutton-Brock et al. dominance index (Clutton-Brock et al. 1979, 1982).
# Remember:
# If A won over B, add 1 to [B,A]
# If A won over C, add 1 ro [C,A]
# -> So the A column represents the species that A dominates (A wins)
# If A lost to B, add 1 to [A,B]
# If A lost to C, add 1 to [A,C]
# -> So the A row represents the species that dominate A (A losses)
# (B + b + 1)/(L + l + 1)
cbi <- function(df, species) {
# number of entries in column 'species' that are neither 0 nor NA:
B <- sum(df[,species] != 0, na.rm = TRUE)
# species corresponding to those entries:
dominates <- rownames(df)[which(df[,species] != 0)]
# number of entries in the columns corresponding to those species that are neither 0 nor NA:
if (length(dominates) > 0) {
b <- sum(subset(df, select = dominates) != 0, na.rm = TRUE)
} else {
b <- 0
}
# number of entries in row B that are neither 0 nor NA:
L <- sum(df[species,], na.rm = TRUE)
# species corresponding to those entries:
dominated <- colnames(df)[which(df[species,] != 0)]
# number of entries in the rows corresponding to those species that are neither 0 nor NA:
if (length(dominated) > 0) {
l <- sum(subset(df, rownames(df) %in% dominated) != 0, na.rm = TRUE)
} else {
l <- 0
}
return((B + b + 1)/(L + l + 1))
}
# Apply it to the full matrix:
cbi_hierarchy <- function(df) {
cbi_matrix <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
colnames(cbi_matrix) <- c("Species", "CBI")
for (i in 1:nrow(df)) {
cbi_matrix$Species[i] <- rownames(df)[i]
cbi_matrix$CBI[i] <- cbi(df, rownames(df)[i])
}
cbi_matrix <- cbi_matrix[order(cbi_matrix$CBI, decreasing = TRUE),]
return(cbi_matrix)
}
cbi_ncdomall <- cbi_hierarchy(nocornet_domall)
cbi_scoredall <- cbi_hierarchy(scoredall)
cbi_ncbothall <- cbi_hierarchy(nocornet_bothall)
cbi_ncooolive <- cbi_hierarchy(nocornet_live)
cbi_ooovideo <- cbi_hierarchy(scoredoooall)
cbi_ncoooboth <- cbi_hierarchy(nocornet_both)
win_freq <- function(df) {
j <- 1
species <- union(levels(df$Species.A), levels(df$Species.B))
freq_table <- data.frame(matrix(NA, nrow = length(species), ncol = 2))
colnames(freq_table) <- c("Species", "WinFreq")
for(i in species) {
if (grepl("*", i, fixed = TRUE) == TRUE) {
next
} else {
freq_table[j,1] <- i
involved <- sum(df$Species.A == i) + sum(df$Species.B == i)
won <- sum(df$Winner == i)
freq_table[j,2] <- won/involved
j <- j + 1
}
}
return(freq_table[complete.cases(freq_table),])
}
Compare the sizes estimated from live observations with those reported by Randall (2005):
# Maximum attained sizes according to Randall (2005):
sizes <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", strip.white = TRUE, header = TRUE)
colnames(sizes)[1] <- "Species"
booksizes <- sizes$Attain[c(1:3,5:8,10:24)]
# Distributions of sizes observed in the field:
fishsize <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FieldFishSize.csv", sep = ",", strip.white = TRUE)
# Calculate the mean, median, and maximum of the observed size distributions:
meansizes <- colMeans(fishsize[,c(2:4,6:9,11:25)], na.rm = TRUE)
colMedians <- function(matrix) {
return(apply(matrix, 2, median, na.rm = TRUE))
}
colMax <- function(matrix) {
return(apply(matrix, 2, max, na.rm = TRUE))
}
mediansizes <- colMedians(fishsize[,c(2:4,6:9,11:25)])
maxsizes <- colMax(fishsize[,c(2:4,6:9,11:25)])
# Do book sizes correlate with live observation sizes?
# a <- cor.test(meansizes, mediansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(a$estimate[[1]],5), ", p-value: ", a$p.value, sep = "")
# b <- cor.test(meansizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(b$estimate[[1]],5), ", p-value: ", b$p.value, sep = "")
# c <- cor.test(mediansizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(c$estimate[[1]],5), ", p-value: ", c$p.value, sep = "")
d <- cor.test(meansizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(d$estimate[[1]],5), ", p-value: ", d$p.value, sep = "")
## [1] "correlation: 0.9969, p-value: 1.4440712484574e-23"
e <- cor.test(meansizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(e$estimate[[1]],5), ", p-value: ", e$p.value, sep = "")
## [1] "correlation: 0.94797, p-value: 2.11724427906556e-11"
f <- cor.test(mediansizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(f$estimate[[1]],5), ", p-value: ", f$p.value, sep = "")
## [1] "correlation: 0.94834, p-value: 1.973173161907e-11"
# Non-unique ranks, gives an error:
# g <- cor.test(booksizes, mediansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(g$estimate[[1]],5), ", p-value: ", g$p.value, sep = "")
# h <- cor.test(booksizes, meansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(h$estimate[[1]],5), ", p-value: ", h$p.value, sep = "")
# i <- cor.test(booksizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(i$estimate[[1]],5), ", p-value: ", i$p.value, sep = "")
j <- cor.test(booksizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(j$estimate[[1]],5), ", p-value: ", j$p.value, sep = "")
## [1] "correlation: 0.83931, p-value: 1.04831251950256e-06"
k <- cor.test(booksizes, meansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(k$estimate[[1]],5), ", p-value: ", k$p.value, sep = "")
## [1] "correlation: 0.85355, p-value: 4.41082818979503e-07"
l <- cor.test(booksizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(l$estimate[[1]],5), ", p-value: ", l$p.value, sep = "")
## [1] "correlation: 0.83313, p-value: 1.48832376900543e-06"
# Reassign the mediansizes vector so that it includes Juvenile Parrotfish:
mediansizes <- colMedians(fishsize[,c(2:4,6:25)])
# Create a new data frame that associates each species with the median of its observed
# size distribution, excluding Cornetfish and "Striated Surgeonfish":
obs_sizes <- data.frame(matrix(ncol = 2, nrow = length(mediansizes)))
colnames(obs_sizes) <- c("Species", "Size")
obs_sizes[,1] <- as.character(sizes$Species[c(1:3,5:24)])
obs_sizes[,2] <- as.numeric(mediansizes)
Include mass (scalled isometrically and allometrically) as an additional variable:
mass <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/body_to_mass_data.csv", sep = ",", strip.white = TRUE)
colnames(mass)[1] <- "Species"
lengthmass <- merge(obs_sizes, mass, by = "Species")
lengthmass[,5] <- lengthmass[,3]*(lengthmass[,2]^lengthmass[,4])
lengthmass[,6] <- lengthmass[,5]^(2/3)
colnames(lengthmass)[c(2,5,6)] <- c("Length", "Mass", "Mass2_3")
Test the relationship between the two measures of dominance (CBI and winning frequency) and body mass or body mass raised to 2/3 by creating a single data frame that contains all the four variables. Note: from this point on, we focus exclusively on the combined (live plus video) dataset of 23 species (no cornetfish) and one-on-one-interactions.
both_test <- merge(cbi_ncoooboth, win_freq(both_ooo4[both_ooo4$Species.A != "Cornetfish" & both_ooo4$Species.B != "Cornetfish",]), by = "Species")
both_test <- merge(both_test, lengthmass, by = "Species")
# Number of interactions:
nrow(both_ooo4[both_ooo4$Species.A != "Cornetfish" & both_ooo4$Species.B != "Cornetfish",])
## [1] 401
cbi_mass <- lm(both_test$CBI ~ both_test$Mass)
summary(cbi_mass)$coefficients[2, 4]
## [1] 0.4754628
summary(cbi_mass)$r.squared
## [1] 0.02452714
winfreq_mass <- lm(both_test$WinFreq ~ both_test$Mass)
summary(winfreq_mass)$coefficients[2, 4]
## [1] 0.1346416
summary(winfreq_mass)$r.squared
## [1] 0.1033787
cbi_mass2_3 <- lm(both_test$CBI ~ both_test$Mass2_3)
summary(cbi_mass2_3)$coefficients[2, 4]
## [1] 0.2865715
summary(cbi_mass2_3)$r.squared
## [1] 0.05387003
winfreq_mass2_3 <- lm(both_test$WinFreq ~ both_test$Mass2_3)
summary(winfreq_mass2_3)$coefficients[2, 4]
## [1] 0.1029628
summary(winfreq_mass2_3)$r.squared
## [1] 0.1215855
# How far is spotted porcupinefish from the median in terms of the interquartile range?
# (Note: this is a simple measure of the extremeness of an outlier)
(both_test$Mass[both_test$Species == "Spotted Porcupinefish"] - median(both_test$Mass))/IQR(both_test$Mass)
## [1] 9.456649
(both_test$Mass2_3[both_test$Species == "Spotted Porcupinefish"] - median(both_test$Mass2_3))/IQR(both_test$Mass2_3)
## [1] 4.972994
# Exclude spotted porcupinefish:
nocornet_nosp_both <- exclude(nocornet_both, "Spotted Porcupinefish")
cbi_ncoooboth_nosp <- cbi_hierarchy(nocornet_nosp_both)
both_test_nosp <- merge(cbi_ncoooboth_nosp, win_freq(both_ooo4[!both_ooo4$Species.A %in% c("Cornetfish", "Spotted Porcupinefish") & !both_ooo4$Species.B %in% c("Cornetfish", "Spotted Porcupinefish"),]), by = "Species")
both_test_nosp <- merge(both_test_nosp, lengthmass, by = "Species")
# Number of data points:
nrow(both_ooo4[!both_ooo4$Species.A %in% c("Cornetfish", "Spotted Porcupinefish") & !both_ooo4$Species.B %in% c("Cornetfish", "Spotted Porcupinefish"),])
## [1] 398
cbi_mass_nosp <- lm(both_test_nosp$CBI ~ both_test_nosp$Mass)
summary(cbi_mass_nosp)$coefficients[2, 4]
## [1] 0.0301074
summary(cbi_mass_nosp)$r.squared
## [1] 0.214146
winfreq_mass_nosp <- lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass)
summary(winfreq_mass_nosp)$coefficients[2, 4]
## [1] 0.2470032
summary(winfreq_mass_nosp)$r.squared
## [1] 0.06638793
cbi_mass2_3_nosp <- lm(both_test_nosp$CBI ~ both_test_nosp$Mass2_3)
summary(cbi_mass2_3_nosp)$coefficients[2, 4]
## [1] 0.03673997
summary(cbi_mass2_3_nosp)$r.squared
## [1] 0.2003125
winfreq_mass2_3_nosp <- lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass2_3)
summary(winfreq_mass2_3_nosp)$coefficients[2, 4]
## [1] 0.2267295
summary(winfreq_mass2_3_nosp)$r.squared
## [1] 0.07216368
This is the code used to generate Figure 1 of the manuscript (NOTE: This figure summarizes the results of the regression analyses uncorrected for phylogenetic relatedness, and is now given as Figure S3 in the Supplementary Information.)
# png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.png", width = 3600, height = 3600, pointsize = 96)
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Supp_Figure_3.png", width = 3600, height = 3600, pointsize = 96)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.eps")
spec_dec <- function(x, k) trimws(format(round(x, k), nsmall = k)) # credit: Jeromy Anglim, Stack Overflow
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass, both_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(lm(both_test$CBI ~ both_test$Mass), lwd = 5)
abline(lm(both_test_nosp$CBI ~ both_test_nosp$Mass), lwd = 5, col = "gray70")
rsquared <- summary(cbi_mass)$r.squared
pvalue <- summary(cbi_mass)$coefficients[2,4]
rsquared_nosp <- summary(cbi_mass_nosp)$r.squared
pvalue_nosp <- summary(cbi_mass_nosp)$coefficients[2,4]
text(800, 10.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(800, 9.1, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(800, 7.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(800, 6, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass2_3, both_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$CBI ~ both_test$Mass2_3), lwd = 5)
abline(lm(both_test_nosp$CBI ~ both_test_nosp$Mass2_3), lwd = 5, col = "gray70")
rsquared <- summary(cbi_mass2_3)$r.squared
pvalue <- summary(cbi_mass2_3)$coefficients[2,4]
rsquared_nosp <- summary(cbi_mass2_3_nosp)$r.squared
pvalue_nosp <- summary(cbi_mass2_3_nosp)$coefficients[2,4]
text(15, 10.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(15, 9.1, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(15, 7.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(15, 6, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass, both_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$WinFreq ~ both_test$Mass), lwd = 5)
abline(lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass), lwd = 5, col = "gray70")
rsquared <- summary(winfreq_mass)$r.squared
pvalue <- summary(winfreq_mass)$coefficients[2,4]
rsquared_nosp <- summary(winfreq_mass_nosp)$r.squared
pvalue_nosp <- summary(winfreq_mass_nosp)$coefficients[2,4]
text(800, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(800, 0.4 - 1.3/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(800, 0.4 - (1.3+1.8)/16, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(800, 0.4 - (2*1.3+1.8)/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass2_3, both_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$WinFreq ~ both_test$Mass2_3), lwd = 5)
abline(lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass2_3), lwd = 5, col = "gray70")
rsquared <- summary(winfreq_mass2_3)$r.squared
pvalue <- summary(winfreq_mass2_3)$coefficients[2,4]
rsquared_nosp <- summary(winfreq_mass2_3_nosp)$r.squared
pvalue_nosp <- summary(winfreq_mass2_3_nosp)$coefficients[2,4]
text(60, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(60, 0.4 - 1.3/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(60, 0.4 - (1.3+1.8)/16, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(60, 0.4 - (2*1.3+1.8)/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()
For each interaction, multiply the body mass of either species (scaled isometrically or allometrically) by its group size or group size squared. See which species maximizes the product, and whether it matches the observed winner. Run a binomial test on the proportion of successful predictions.
both_group <- both4[both4$Group.A != both4$Group.B,]
# Clean-up: drop rows containing NAs:
both_group <- both_group[complete.cases(both_group),]
# Exclude cornetfish:
both_group <- both_group[both_group$Species.A != "Cornetfish" &
both_group$Species.B != "Cornetfish",]
# Number of data points:
nrow(both_group)
## [1] 83
# Get the body sizes for both competitors:
for(i in 1:nrow(both_group)) {
species_a <- as.character(both_group$Species.A[i])
species_b <- as.character(both_group$Species.B[i])
both_group[i,6] <- lengthmass$Mass[lengthmass$Species == species_a]
both_group[i,7] <- lengthmass$Mass[lengthmass$Species == species_b]
}
colnames(both_group)[6:7] <- c("BodySize.A", "BodySize.B")
# Convert group sizes to the right object type (from 'character' to 'double'):
both_group$Group.A <- as.numeric(both_group$Group.A)
both_group$Group.B <- as.numeric(both_group$Group.B)
# Make predictions based on the linear law:
for(i in 1:nrow(both_group)) {
if (both_group$Group.A[i]*both_group$BodySize.A[i] > both_group$Group.B[i]*both_group$BodySize.B[i]) {
both_group[i,8] <- both_group$Species.A[i]
} else {
both_group[i,8] <- both_group$Species.B[i]
}
}
# Make predictions based on the square law:
for(i in 1:nrow(both_group)) {
if ((both_group$Group.A[i]^2)*both_group$BodySize.A[i] > (both_group$Group.B[i]^2)*both_group$BodySize.B[i]) {
both_group[i,9] <- both_group$Species.A[i]
} else {
both_group[i,9] <- both_group$Species.B[i]
}
}
colnames(both_group)[8:9] <- c("Linear.Law.Prediction", "Square.Law.Prediction")
# Test if the linear law works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Linear.Law.Prediction[i])) == TRUE) {
both_group[i,10] <- 1
} else {
both_group[i,10] <- 0
}
}
# Test if the square law works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Square.Law.Prediction[i])) == TRUE) {
both_group[i,11] <- 1
} else {
both_group[i,11] <- 0
}
}
colnames(both_group)[10:11] <- c("Linear.True", "Square.True")
# Proportion of successful predictions:
sum(both_group$Linear.True)/nrow(both_group)
## [1] 0.6626506
sum(both_group$Square.True)/nrow(both_group)
## [1] 0.4939759
# Test significance:
dbinom(sum(both_group$Linear.True), size = nrow(both_group), prob = 1/2)
## [1] 0.001053886
dbinom(sum(both_group$Square.True), size = nrow(both_group), prob = 1/2)
## [1] 0.08679764
# When do the two laws give different predictions?
dif <- both_group[both_group$Linear.Law.Prediction != both_group$Square.Law.Prediction,]
dbinom(sum(dif$Linear.True), size = nrow(dif), prob = 1/2)
## [1] 0.001087189
# Now repeat all of this for body mass raised to two thirds:
both_group2 <- both4[both4$Group.A != both4$Group.B,]
both_group2 <- both_group2[complete.cases(both_group2),]
both_group2 <- both_group2[both_group2$Species.A != "Cornetfish" &
both_group2$Species.B != "Cornetfish",]
for(i in 1:nrow(both_group2)) {
species_a <- as.character(both_group2$Species.A[i])
species_b <- as.character(both_group2$Species.B[i])
both_group2[i,6] <- lengthmass$Mass2_3[lengthmass$Species == species_a]
both_group2[i,7] <- lengthmass$Mass2_3[lengthmass$Species == species_b]
}
colnames(both_group2)[6:7] <- c("BodySize.A", "BodySize.B")
both_group2$Group.A <- as.numeric(both_group2$Group.A)
both_group2$Group.B <- as.numeric(both_group2$Group.B)
for(i in 1:nrow(both_group2)) {
if (both_group2$Group.A[i]*both_group2$BodySize.A[i] > both_group2$Group.B[i]*both_group2$BodySize.B[i]) {
both_group2[i,8] <- both_group2$Species.A[i]
} else {
both_group2[i,8] <- both_group2$Species.B[i]
}
}
for(i in 1:nrow(both_group2)) {
if ((both_group2$Group.A[i]^2)*both_group2$BodySize.A[i] > (both_group2$Group.B[i]^2)*both_group2$BodySize.B[i]) {
both_group2[i,9] <- both_group2$Species.A[i]
} else {
both_group2[i,9] <- both_group2$Species.B[i]
}
}
colnames(both_group2)[8:9] <- c("Linear.Law.Prediction", "Square.Law.Prediction")
for(i in 1:nrow(both_group2)) {
if (identical(as.character(both_group2$Winner[i]), as.character(both_group2$Linear.Law.Prediction[i])) == TRUE) {
both_group2[i,10] <- 1
} else {
both_group2[i,10] <- 0
}
}
for(i in 1:nrow(both_group2)) {
if (identical(as.character(both_group2$Winner[i]), as.character(both_group2$Square.Law.Prediction[i])) == TRUE) {
both_group2[i,11] <- 1
} else {
both_group2[i,11] <- 0
}
}
colnames(both_group2)[10:11] <- c("Linear.True", "Square.True")
sum(both_group2$Linear.True)/nrow(both_group2)
## [1] 0.5662651
sum(both_group2$Square.True)/nrow(both_group2)
## [1] 0.3614458
dbinom(sum(both_group2$Linear.True), size = nrow(both_group2), prob = 1/2)
## [1] 0.04240454
dbinom(sum(both_group2$Square.True), size = nrow(both_group2), prob = 1/2)
## [1] 0.003597748
dif2 <- both_group2[both_group2$Linear.Law.Prediction != both_group2$Square.Law.Prediction,]
dbinom(sum(dif2$Linear.True), size = nrow(dif2), prob = 1/2)
## [1] 0.0001001358
Plot a curve showing how many unique interactions you get with each new video scored:
curve_analysis <- drop_ufos(scored)
curve_analysis <- curve_analysis[complete.cases(curve_analysis),]
curve_analysis <- curve_analysis[,c(4,5,6,10,13)]
# Convert the dates and times into a readable format and sort the dataset by them:
library(chron)
curve_analysis$Date.of.recording <- as.Date(curve_analysis$Date.of.recording, format = "%d %B %y")
# This will makes the times readable by 'chron':
curve_analysis$Time.of.Recording <- paste(curve_analysis$Time.of.Recording, ":00", sep = "")
curve_analysis$Time.of.Recording <- chron(times = curve_analysis$Time.of.Recording)
# Sort by date and time:
curve_analysis <- curve_analysis[order(curve_analysis$Date.of.recording, curve_analysis$Time.of.Recording, decreasing = TRUE),]
species_pairs <- list()
rarefaction <- data.frame(matrix(NA, nrow = nrow(curve_analysis), ncol = 2))
colnames(rarefaction) <- c("Observation", "Unique Pairs")
for(i in 1:nrow(curve_analysis)) {
new_pair <- c(as.character(curve_analysis$Species.A[i]), as.character(curve_analysis$Species.B[i]))
species_pairs <- c(species_pairs, list(new_pair))
rarefaction[i,1] <- i
rarefaction[i,2] <- length(unique(species_pairs))
}
plot(rarefaction$Observation, rarefaction$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species")
# Now order the observations according to scoring time rather than recording time:
curve_analysis2 <- drop_ufos(scored)
curve_analysis2 <- curve_analysis2[complete.cases(curve_analysis2),]
curve_analysis2 <- curve_analysis2[,c(3,10,13)]
curve_analysis2$Date.of.scoring <- as.Date(curve_analysis2$Date.of.scoring, format = "%d %B %Y")
curve_analysis2 <- curve_analysis2[order(curve_analysis2$Date.of.scoring, decreasing = TRUE),]
species_pairs2 <- list()
rarefaction2 <- data.frame(matrix(NA, nrow = nrow(curve_analysis), ncol = 2))
colnames(rarefaction2) <- c("Observation", "Unique Pairs")
for(i in 1:nrow(curve_analysis2)) {
new_pair <- c(as.character(curve_analysis2$Species.A[i]), as.character(curve_analysis2$Species.B[i]))
species_pairs2 <- c(species_pairs2, list(new_pair))
rarefaction2[i,1] <- i
rarefaction2[i,2] <- length(unique(species_pairs2))
}
plot(rarefaction2$Observation, rarefaction2$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species", main = "Contribution of new observations to dataset expansion")
This is the code used to generate Figure S4:
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_S1.png", width = 2880, height = 1800, pointsize = 64)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_S1.eps")
plot(rarefaction2$Observation, rarefaction2$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species", main = "Contribution of new observations to dataset expansion")
dev.off()
diet <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/diet.csv", sep = ",", strip.white = TRUE)
# Add dietary information for both Species A and Species B
for(i in 1:nrow(both_group)) {
both_group[i, 12] <- diet$Diet[diet$Species == as.character(both_group$Species.A[i])]
both_group[i, 13] <- diet$Diet[diet$Species == as.character(both_group$Species.B[i])]
}
colnames(both_group)[12:13] <- c("Diet.A", "Diet.B")
# Now subsample the test frame to include only those interactions that take place between
# species with the same diet:
test_benth <- both_group[both_group$Diet.A == "Benthic Invertebrate Consumer" &
both_group$Diet.B == "Benthic Invertebrate Consumer",]
nrow(test_benth)
## [1] 0
test_omni <- both_group[both_group$Diet.A == "Omnivore" & both_group$Diet.B == "Omnivore",]
nrow(test_omni)
## [1] 16
test_herb <- both_group[both_group$Diet.A == "Herbivore/Detritivore" & both_group$Diet.B == "Herbivore/Detritivore",]
nrow(test_herb)
## [1] 4
test_coral <- both_group[both_group$Diet.A == "Corallivore" & both_group$Diet.B == "Corallivore",]
nrow(test_coral)
## [1] 0
# Test the validity of Lanchester's laws for interaction subsets of non-zero length.
# Omnivores:
sum(test_omni$Linear.True)/nrow(test_omni)
## [1] 0.8125
sum(test_omni$Square.True)/nrow(test_omni)
## [1] 0.5625
dbinom(sum(test_omni$Linear.True), size = nrow(test_omni), prob = 1/2)
## [1] 0.008544922
dbinom(sum(test_omni$Square.True), size = nrow(test_omni), prob = 1/2)
## [1] 0.1745605
dif3 <- test_omni[test_omni$Linear.Law.Prediction != test_omni$Square.Law.Prediction,]
nrow(dif3)
## [1] 4
sum(dif3$Linear.True)/nrow(dif3)
## [1] 1
dbinom(sum(dif3$Linear.True), size = nrow(dif3), prob = 1/2)
## [1] 0.0625
# Herbivores/Detritivores:
sum(test_herb$Linear.True)/nrow(test_herb)
## [1] 1
sum(test_herb$Square.True)/nrow(test_herb)
## [1] 0.25
dbinom(sum(test_herb$Linear.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.0625
dbinom(sum(test_herb$Square.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.25
dif4 <- test_herb[test_herb$Linear.Law.Prediction != test_herb$Square.Law.Prediction,]
nrow(dif4)
## [1] 3
sum(dif4$Linear.True)/nrow(dif4)
## [1] 1
dbinom(sum(dif4$Linear.True), size = nrow(dif4), prob = 1/2)
## [1] 0.125
# Re-run everything with mass raised to two thirds:
for(i in 1:nrow(both_group2)) {
both_group2[i, 12] <- diet$Diet[diet$Species == as.character(both_group2$Species.A[i])]
both_group2[i, 13] <- diet$Diet[diet$Species == as.character(both_group2$Species.B[i])]
}
colnames(both_group2)[12:13] <- c("Diet.A", "Diet.B")
test_benth2 <- both_group2[both_group2$Diet.A == "Benthic Invertebrate Consumer" &
both_group2$Diet.B == "Benthic Invertebrate Consumer",]
test_omni2 <- both_group2[both_group2$Diet.A == "Omnivore" & both_group2$Diet.B == "Omnivore",]
test_herb2 <- both_group2[both_group2$Diet.A == "Herbivore/Detritivore" & both_group2$Diet.B == "Herbivore/Detritivore",]
sum(test_omni2$Linear.True)/nrow(test_omni2)
## [1] 0.5625
sum(test_omni2$Square.True)/nrow(test_omni2)
## [1] 0.25
dbinom(sum(test_omni2$Linear.True), size = nrow(test_omni2), prob = 1/2)
## [1] 0.1745605
dbinom(sum(test_omni2$Square.True), size = nrow(test_omni2), prob = 1/2)
## [1] 0.027771
dif5 <- test_omni2[test_omni2$Linear.Law.Prediction != test_omni2$Square.Law.Prediction,]
nrow(dif5)
## [1] 5
sum(dif5$Linear.True)/nrow(dif5)
## [1] 1
dbinom(sum(dif5$Linear.True), size = nrow(dif5), prob = 1/2)
## [1] 0.03125
sum(test_herb2$Linear.True)/nrow(test_herb2)
## [1] 0.75
sum(test_herb2$Square.True)/nrow(test_herb2)
## [1] 0
dbinom(sum(test_herb2$Linear.True), size = nrow(test_herb2), prob = 1/2)
## [1] 0.25
dbinom(sum(test_herb2$Square.True), size = nrow(test_herb2), prob = 1/2)
## [1] 0.0625
dif6 <- test_herb2[test_herb2$Linear.Law.Prediction != test_herb2$Square.Law.Prediction,]
nrow(dif6)
## [1] 3
sum(dif6$Linear.True)/nrow(dif6)
## [1] 1
dbinom(sum(dif6$Linear.True), size = nrow(dif6), prob = 1/2)
## [1] 0.125
These include exploratory analyses as well as applications of previously proposed methods to the fish dataset (e.g., Shelley et al.’s (2004) logistic regression).
Subsample the data sheets down to the 10 most often interacting species. This will increase completeness (substantially) and linearity (slightly), but give you fewer data points to work with.
mostoften <- function(sheet, df, n) {
occurrences <- tail(sort(table(unlist(sheet[,c("Species.A", "Species.B")]))), n)
max_n <- as.character(as.data.frame(occurrences)[,1])
max_n_matrix <- subset(df, select = max_n)
max_n_matrix <- subset(max_n_matrix, rownames(max_n_matrix) %in% max_n)
max_n_matrix <- max_n_matrix[order(row.names(max_n_matrix)), order(colnames(max_n_matrix))]
return(max_n_matrix)
}
dommatrix_max10 <- mostoften(sheet4, dommatrix, 10)
scoredmatrix_max10 <- mostoften(scored4, scoredmatrix, 10)
bothmatrix_max10 <- mostoften(both4, bothmatrix, 10)
domooo_max10 <- mostoften(sheet_ooo4, domooo, 10)
scoredooo_max10 <- mostoften(scored_ooo4, scoredooo, 10)
bothooo_max10 <- mostoften(both_ooo4, bothooo, 10)
dommax10 <- convert_to_df(dommatrix_max10)
scoredmax10 <- convert_to_df(scoredmatrix_max10)
bothmax10 <- convert_to_df(bothmatrix_max10)
liveooomax10 <- convert_to_df(domooo_max10)
scoredooomax10 <- convert_to_df(scoredooo_max10)
bothooomax10 <- convert_to_df(bothooo_max10)
paste("The completeness of the live-only 10-species matrix is ", compl(dommax10), "%", sep = "")
paste("The completeness of the video-only 10-species matrix is ", compl(scoredmax10), "%", sep = "")
paste("The completeness of the combined 10-species matrix is ", compl(bothmax10), "%", sep = "")
paste("The completeness of the live-only one-on-one 10-species matrix is ", compl(liveooomax10), "%", sep = "")
paste("The completeness of the video-only one-on-one 10-species matrix is ", compl(scoredooomax10), "%", sep = "")
paste("The completeness of the combined one-on-one 10-species matrix is ", compl(bothooomax10), "%", sep = "")
cbi_max10 <- cbi_hierarchy(dommax10)
cbi_scoredmax10 <- cbi_hierarchy(scoredmax10)
cbi_bothmax10 <- cbi_hierarchy(bothmax10)
cbi_liveooomax10 <- cbi_hierarchy(liveooomax10)
cbi_scoredooomax10 <- cbi_hierarchy(scoredooomax10)
cbi_bothooomax10 <- cbi_hierarchy(bothooomax10)
devries(dommax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(scoredmax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(bothmax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(liveooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(scoredooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(bothooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
A third measure of dominance that is particularly applicable to weakly but significantly linear dominance hierarchies, implemented as described in Schmid & de Vries (2013).
library(compete)
i_and_si_liveooo <- isi13(liveoooall, nTries = 10)
rev(i_and_si_liveooo$best_order) # The order is given from least to most dominant by default
i_and_si_videoooo <- isi13(scoredoooall, nTries = 10)
rev(i_and_si_videoooo$best_order)
i_and_si_bothooo <- isi13(bothoooall, nTries = 10)
rev(i_and_si_bothooo$best_order)
i_and_si <- isi13(dommax10, nTries = 10)
rev(i_and_si$best_order)
i_and_si_scored <- isi13(scoredmax10, nTries = 10)
rev(i_and_si_scored$best_order)
i_and_si_both <- isi13(bothmax10, nTries = 10)
rev(i_and_si_both$best_order)
# Create a data frame with alphabetically ordered species and their corresponding I&SI ranks:
isi_both <- i_and_si_bothooo$best_order
isi_both <- setdiff(isi_both, "Cornetfish")
test_isi <- data.frame(matrix(NA, nrow = length(isi_both), ncol = 2))
colnames(test_isi) <- c("Species", "I_SI_Rank")
test_isi$Species <- sort(isi_both)
test_isi$I_SI_Rank <- match(sort(isi_both), isi_both)
Do simple nonlinear models (exponential, logarithmic, quadratic) explain the relationship between individual fighting ability (as determined by CBI and the frequency of winning) and body mass?
# Add mass squared to the test frame:
both_test[,9] <- both_test$Mass^2
colnames(both_test)[9] <- "Mass_Squared"
stats <- function(df) {
results <- data.frame(matrix(NA, nrow = 2, ncol = 14))
# Columns represent different statistics associated with each model. Spearman's rho
# and Pearson's coefficient are included for comparison; if the former is significant
# but the latter is not, this shows that the relationship is monotonic but nonlinear.
colnames(results) <- c("rho", "Spear_p-value", "cor", "Pear_p-value", "linAIC", "expRsq", "exp_p-value", "expAIC", "logRsq", "log_p-value", "logAIC", "quadrRsq", "quadr_p-value", "quadrAIC")
# Rows represent dependent variables:
row.names(results) <- c("CBI_Mass", "WinFreq_Mass")
spear_mass_CBI <- cor.test(df$CBI, df$Mass, method = "spearman", alternative = "two.sided")
spear_mass_WinFreq <- cor.test(df$WinFreq, df$Mass, method = "spearman", alternative = "two.sided")
pear_mass_CBI <- cor.test(df$CBI, df$Mass, method = "pearson", alternative = "two.sided")
pear_mass_WinFreq <- cor.test(df$WinFreq, df$Mass, method = "pearson", alternative = "two.sided")
linear_mass_CBI <- lm(df$CBI ~ df$Mass)
expo_mass_CBI <- lm(log(df$CBI) ~ df$Mass)
loga_mass_CBI <- lm(df$CBI ~ log(df$Mass))
quadra_mass_CBI <- lm(df$CBI ~ df$Mass + df$Mass_Squared)
linear_mass_WinFreq <- lm(df$WinFreq ~ df$Mass)
# No exponential model for winning frequency: the minimum winning frequency observed
# in the dataset is 0, and simple log transformation therefore cannot be applied to it.
loga_mass_WinFreq <- lm(df$WinFreq ~ log(df$Mass))
quadra_mass_WinFreq <- lm(df$WinFreq ~ df$Mass + df$Mass_Squared)
results[1,1] <- spear_mass_CBI$estimate; results[1,2] <- spear_mass_CBI$p.value
results[1,3] <- pear_mass_CBI$estimate; results[1,4] <- pear_mass_CBI$p.value
results[1,5] <- AIC(linear_mass_CBI); results[1,6] <- summary(expo_mass_CBI)$r.squared
results[1,7] <- summary(expo_mass_CBI)$coefficients[,4][[2]]
results[1,8] <- AIC(expo_mass_CBI); results[1,9] <- summary(loga_mass_CBI)$r.squared
results[1,10] <- summary(loga_mass_CBI)$coefficients[,4][[2]]
results[1,11] <- AIC(loga_mass_CBI); results[1,12] <- summary(quadra_mass_CBI)$r.squared
results[1,13] <- summary(quadra_mass_CBI)$coefficients[,4][[2]]
results[1,14] <- AIC(quadra_mass_CBI)
results[2,1] <- spear_mass_WinFreq$estimate; results[2,2] <- spear_mass_WinFreq$p.value
results[2,3] <- pear_mass_WinFreq$estimate; results[2,4] <- pear_mass_WinFreq$p.value
results[2,5] <- AIC(linear_mass_WinFreq)
results[2,9] <- summary(loga_mass_WinFreq)$r.squared
results[2,10] <- summary(loga_mass_WinFreq)$coefficients[,4][[2]]
results[2,11] <- AIC(loga_mass_WinFreq)
results[2,12] <- summary(quadra_mass_WinFreq)$r.squared
results[2,13] <- summary(quadra_mass_WinFreq)$coefficients[,4][[2]]
results[2,14] <- AIC(quadra_mass_WinFreq)
return(results)
}
both_test_results <- stats(both_test)
# Adjust the p-values for multiple comparisons:
library(stats)
pvalues <- unlist(both_test_results[,c(2,4,7,10,13)])
p.adjust(pvalues, "bonferroni", sum(!is.na(pvalues)))
The code below is an implementation of the logistic regression analysis described by Shelley et al. (2004). Given the paucity of data, no pairwise comparisons were made; each of the selected species was compared against all other species it interacted with. Unlike Shelley et al. (2004), we regress against both group size and group size squared.
# Extract those species that participated in an interaction with more than one individual
# at least a given number of times:
atleastoncemorethan <- function(df, n) {
output <- vector()
species <- union(unique(df$Species.A), unique(df$Species.B))
for (i in species) {
a.frame <- df[df$Species.A == i,]
b.frame <- df[df$Species.B == i,]
if (sum(!is.na(a.frame$Group.A) & a.frame$Group.A > 1) + sum(!is.na(b.frame$Group.B) & b.frame$Group.B > 1) >= n) {
output <- c(output, as.character(i))
}
}
return(unique(output))
}
# Some ad hoc data clean-up before applying the atleastoncemorethan() function:
scored4$Group.B <- as.numeric(as.character(scored4$Group.B))
both4$Group.B <- as.numeric(both4$Group.B)
# Include only those species that had a group size greater than 1 at least twice. This
# ensures that there is variability in the predictor variable.
live_logres <- sheet4[sheet4$Species.A %in% atleastoncemorethan(sheet4, 2) &
sheet4$Species.B %in% atleastoncemorethan(sheet4, 2),]
video_logres <- scored4[scored4$Species.A %in% atleastoncemorethan(scored4, 2) &
scored4$Species.B %in% atleastoncemorethan(scored4, 2),]
both_logres <- both4[both4$Species.A %in% atleastoncemorethan(both4, 2) &
both4$Species.B %in% atleastoncemorethan(both4, 2),]
# Convert from wide to long format: every two successive rows correspond to one
# interaction, with the winner coded as "1" and the loser coded as "0".
wide_to_long <- function(df) {
longframe <- data.frame(matrix(NA, ncol = 3, nrow = 2*nrow(df)))
colnames(longframe) <- c("Species", "Group_Size", "Is_Winner")
for (i in 1:nrow(df)) {
longframe$Species[2*i-1] <- as.character(df$Species.A[i])
longframe$Species[2*i] <- as.character(df$Species.B[i])
longframe$Group_Size[2*i-1] <- df$Group.A[i]
longframe$Group_Size[2*i] <- df$Group.B[i]
if (identical(as.character(df$Species.A[i]), as.character(df$Winner[i])) == TRUE) {
longframe$Is_Winner[2*i-1] <- 1
longframe$Is_Winner[2*i] <- 0
}
if (identical(as.character(df$Species.B[i]), as.character(df$Winner[i])) == TRUE) {
longframe$Is_Winner[2*i-1] <- 0
longframe$Is_Winner[2*i] <- 1
}
}
return(longframe)
}
# Make sure winner is always either A or B, and that you have group sizes for all rows:
live_long <- wide_to_long(live_logres)
nrow(live_long) == nrow(live_long[!is.na(live_long$Is_Winner),])
nrow(live_long) == nrow(live_long[!is.na(live_long$Group_Size),])
video_long <- wide_to_long(video_logres)
nrow(video_long) == nrow(video_long[!is.na(video_long$Is_Winner),])
nrow(video_long) == nrow(video_long[!is.na(video_long$Group_Size),])
both_long <- wide_to_long(both_logres)
nrow(both_long) == nrow(both_long[!is.na(both_long$Is_Winner),])
nrow(both_long) == nrow(both_long[!is.na(both_long$Group_Size),])
# Add group size squared
square_it <- function(df) {
a <- ncol(df)
df[,a+1] <- df$Group_Size^2
colnames(df)[a+1] <- "Group_Size_Squared"
return(df)
}
live_long <- square_it(live_long)
video_long <- square_it(video_long)
both_long <- square_it(both_long)
# Extract those species that won and lost at least 1 interaction. This ensures that there
# is variability in the dependent variable:
won_once <- function(df) {
output <- vector()
species <- unique(df$Species)
for(i in species) {
subset.frame <- df[df$Species == i,]
if(0 %in% subset.frame$Is_Winner & 1 %in% subset.frame$Is_Winner) {
output <- c(output, i)
}
}
return(df[df$Species %in% output,])
}
# This function creates a data frame with logistic regression results. The regression is
# run against group size if x is FALSE and against group size squared if x is TRUE.
library(pscl)
species_logistic <- function(df, x = TRUE) {
j <- 1
species <- unique(df$Species)
results <- data.frame(matrix(NA, nrow = length(species), ncol = 6))
colnames(results) <- c("Species", "Coefficient", "p-value", "McFaddenpseudoRsq", "AIC", "n")
for(i in species) {
if (x == FALSE) {
model <- glm(df$Is_Winner[df$Species == i] ~ df$Group_Size[df$Species == i])
} else {
model <- glm(df$Is_Winner[df$Species == i] ~ df$Group_Size_Squared[df$Species == i] + df$Group_Size[df$Species == i])
}
results[j,1] <- i
results[j,2] <- model$coefficients[2]
results[j,3] <- coef(summary(model))[,4][2]
results[j,4] <- pR2(model)[4]
results[j,5] <- AIC(model)
results[j,6] <- nrow(df[df$Species == i,])
j <- j + 1
}
return(results)
}
# Linear group size
species_logistic(won_once(live_long), FALSE)
species_logistic(won_once(video_long), FALSE)
species_logistic(won_once(both_long), FALSE)
# Group size squared
species_logistic(won_once(live_long), TRUE)
species_logistic(won_once(video_long), TRUE)
species_logistic(won_once(both_long), TRUE)
This is an alternative way of testing Lanchester’s laws. Instead of applying the binomial test to a binary variable (prediction: correct, incorrect), we apply \(\chi^2\) to a 4-entry contingency table (correctly predicted win, incorrectly predicted win, correctly predicted loss, incorrectly predicted loss):
populate_cont_table <- function(df, law) {
# Create a contingency table whose rows represent (1) Wins, (2) Losses predicted
# by the linear law and whose columns represent observed (1) Wins and (2) Losses
cont_table <- matrix(0, nrow = 2, ncol = 2)
# Focus solely on Species A. This is arbitrary; we could instead focus on B by
# interchanging individual rows and individual columns.
for(i in 1:nrow(df)) {
# A won:
if (identical(as.character(df$Species.A[i]), as.character(df$Winner[i])) == TRUE) {
# A was predicted to win:
if (identical(as.character(df$Species.A[i]), as.character(df[i, paste(law, ".Law.Prediction", sep = "")])) == TRUE) {
inc(cont_table[1,1])
# A was predicted to lose:
} else {
inc(cont_table[2,1])
}
# A lost:
} else {
# A was predicted to win:
if (identical(as.character(df$Species.A[i]), as.character(df[i, paste(law, ".Law.Prediction", sep = "")])) == TRUE) {
inc(cont_table[1,2])
# A was predicted to lose:
} else {
inc(cont_table[2,2])
}
}
}
return(cont_table)
}
chisq.test(populate_cont_table(both_group, "Linear"))
chisq.test(populate_cont_table(both_group, "Square"))
chisq.test(populate_cont_table(both_group2, "Linear"))
chisq.test(populate_cont_table(both_group2, "Square"))
# To distinguish between the linear and square laws, we run another chi-squared test.
# This time, the contingency table has the following structure: [1,1] correctly
# predicted by both models, [1,2] correctly predicted by the linear law but not the
# square law, [2,1] correctly predicted by the square law but not the linear law,
# [2,2] predicted incorrectly by both laws.
law_vs_law <- function(df) {
cont_table <- matrix(0, nrow = 2, ncol = 2)
for(i in 1:nrow(df)) {
if (df$Linear.True[i] == 1 & df$Square.True[i] == 1) {
inc(cont_table[1,1])
} else if (df$Linear.True[i] == 1 & df$Square.True[i] == 0) {
inc(cont_table[1,2])
} else if (df$Linear.True[i] == 0 & df$Square.True[i] == 1) {
inc(cont_table[2,1])
} else {
inc(cont_table[2,2])
}
}
return(cont_table)
}
laws_mass <- law_vs_law(both_group)
laws_mass2_3 <- law_vs_law(both_group2)
chisq.test(laws_mass)
chisq.test(laws_mass2_3)
We discovered a discrepancy between the identification of the longnose butterflyfish in our fish size data (FishSize.csv) and Table 1 of our manuscript: the same species was identified as Forcipiger longirostris in the former and as the closely related F. flavissimus in the latter. Based on the photograph we took of the species, we concluded that the latter identification was correct, as the lower half of its eye was silver rather than black (http://reefs.com/2015/12/31/long-nosed-butterflies-part-2-forcipiger/). We therefore replaced the maximum attained size reported for F. longirostris by Nelson (2005) with that reported for F. flavissimus. Below, we recalculate the correlation between the reported attained sizes and the observed sizes:
# Maximum attained sizes according to Randall (2005):
sizes <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", header = TRUE)
colnames(sizes)[1] <- "Species"
booksizes <- sizes$Attain[c(1:3,5:8,10:24)]
# Distributions of sizes observed in the field:
fishsize <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FieldFishSize.csv", sep = ",", strip.white = TRUE)
# Calculate the mean, median, and maximum of the observed size distributions:
meansizes <- colMeans(fishsize[,c(2:4,6:9,11:25)], na.rm = TRUE)
colMedians <- function(matrix) {
return(apply(matrix, 2, median, na.rm = TRUE))
}
colMax <- function(matrix) {
return(apply(matrix, 2, max, na.rm = TRUE))
}
mediansizes <- colMedians(fishsize[,c(2:4,6:9,11:25)])
maxsizes <- colMax(fishsize[,c(2:4,6:9,11:25)])
# Do book sizes correlate with live observation sizes?
g <- cor.test(booksizes, mediansizes, method = "spearman", alternative = "two.sided")
## Warning in cor.test.default(booksizes, mediansizes, method = "spearman", :
## Cannot compute exact p-value with ties
paste("correlation: ", signif(g$estimate[[1]],5), ", p-value: ", g$p.value, sep = "")
## [1] "correlation: 0.71489, p-value: 0.000184888348145873"
h <- cor.test(booksizes, meansizes, method = "spearman", alternative = "two.sided")
## Warning in cor.test.default(booksizes, meansizes, method = "spearman",
## alternative = "two.sided"): Cannot compute exact p-value with ties
paste("correlation: ", signif(h$estimate[[1]],5), ", p-value: ", h$p.value, sep = "")
## [1] "correlation: 0.76279, p-value: 3.66099985642471e-05"
i <- cor.test(booksizes, maxsizes, method = "spearman", alternative = "two.sided")
## Warning in cor.test.default(booksizes, maxsizes, method = "spearman",
## alternative = "two.sided"): Cannot compute exact p-value with ties
paste("correlation: ", signif(i$estimate[[1]],5), ", p-value: ", i$p.value, sep = "")
## [1] "correlation: 0.81934, p-value: 3.09742996782284e-06"
j <- cor.test(booksizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(j$estimate[[1]],5), ", p-value: ", j$p.value, sep = "")
## [1] "correlation: 0.83931, p-value: 1.04831251950256e-06"
k <- cor.test(booksizes, meansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(k$estimate[[1]],5), ", p-value: ", k$p.value, sep = "")
## [1] "correlation: 0.85355, p-value: 4.41082818979503e-07"
l <- cor.test(booksizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(l$estimate[[1]],5), ", p-value: ", l$p.value, sep = "")
## [1] "correlation: 0.83313, p-value: 1.48832376900543e-06"
Going back to Clutton-Brock et al. (1979), we discovered that our original implementation of the CBI-calculating function had an error. When counting the species beaten by those species that were dominated by the subject, we need to exclude the subject itself; the same is true when tallying the species that dominated those species which dominated the subject. This was not done in the original code, but it is done below:
# Remember:
# The A column represents the species that A dominates (A wins)
# The A row represents the species that dominate A (A losses)
# (B + b + 1)/(L + l + 1)
cbi_new <- function(df, species) {
# number of entries in column 'species' that are neither 0 nor NA:
B <- sum(df[,species] != 0, na.rm = TRUE)
# species corresponding to those entries:
dominates <- rownames(df)[which(df[,species] != 0)]
# number of entries in the columns corresponding to those species that are neither 0 nor NA:
if (length(dominates) > 0) {
# exclude 'species' itself from consideration
b <- sum(df[rownames(df) != species, colnames(df) %in% species] != 0, na.rm = TRUE)
} else {
b <- 0
}
# number of entries in row B that are neither 0 nor NA:
L <- sum(df[species,], na.rm = TRUE)
# species corresponding to those entries:
dominated <- colnames(df)[which(df[species,] != 0)]
# number of entries in the rows corresponding to those species that are neither 0 nor NA:
if (length(dominated) > 0) {
# exclude 'species' itself from consideration
l <- sum(df[rownames(df) %in% dominated, colnames(df) != species] != 0, na.rm = TRUE)
} else {
l <- 0
}
return((B + b + 1)/(L + l + 1))
}
# Apply it to the full matrix:
cbi_hierarchy_new <- function(df) {
cbi_matrix <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
colnames(cbi_matrix) <- c("Species", "CBI")
for (i in 1:nrow(df)) {
cbi_matrix$Species[i] <- rownames(df)[i]
cbi_matrix$CBI[i] <- cbi_new(df, rownames(df)[i])
}
cbi_matrix <- cbi_matrix[order(cbi_matrix$CBI, decreasing = TRUE),]
return(cbi_matrix)
}
First, we need to drop the “Juvenile Parrotfish” category from the analysis, since it is a mixture of multiple species:
nojp_domall <- exclude(nocornet_domall, "Juvenile Parrotfish")
nojp_scoredall <- exclude(scoredall, "Juvenile Parrotfish")
nojp_bothall <- exclude(nocornet_bothall, "Juvenile Parrotfish")
nojp_live <- exclude(nocornet_live, "Juvenile Parrotfish")
nojp_scored <- exclude(scoredoooall, "Juvenile Parrotfish")
nojp_both <- exclude(nocornet_both, "Juvenile Parrotfish")
# Recalculate CBI:
cbi_nojp_domall <- cbi_hierarchy_new(nojp_domall)
cbi_nojp_scoredall <- cbi_hierarchy_new(nojp_scoredall)
cbi_nojp_bothall <- cbi_hierarchy_new(nojp_bothall)
cbi_nojp_ooolive <- cbi_hierarchy_new(nojp_live)
cbi_nojp_ooovideo <- cbi_hierarchy_new(nojp_scored)
cbi_nojp_oooboth <- cbi_hierarchy_new(nojp_both)
# Recalculate the frequency of winning and combine the results with the CBI hierarchy
# and mass-length data into a single data frame for further testing. We will again focus
# on one-on-one data only from now on.
# This is a generalized version of the no_surgeon() function defined above:
drop_taxa <- function(df, taxa_to_drop) {
return(df[!df$Species.A %in% taxa_to_drop & !df$Species.B %in% taxa_to_drop,])
}
sheet_ooo5 <- drop_taxa(sheet_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
scored_ooo5 <- drop_taxa(scored_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
both_ooo5 <- drop_taxa(both_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
sheet_phylo_test <- merge(cbi_nojp_ooolive, win_freq(sheet_ooo5), by = "Species")
sheet_phylo_test <- merge(sheet_phylo_test, lengthmass, by = "Species")
# We need to provide row names to ensure that each row of the data frame is linked to the
# correct tree tip:
rownames(sheet_phylo_test) <- sheet_phylo_test$Species
scored_phylo_test <- merge(cbi_nojp_ooovideo, win_freq(scored_ooo5), by = "Species")
scored_phylo_test <- merge(scored_phylo_test, lengthmass, by = "Species")
rownames(scored_phylo_test) <- scored_phylo_test$Species
both_phylo_test <- merge(cbi_nojp_oooboth, win_freq(both_ooo5), by = "Species")
both_phylo_test <- merge(both_phylo_test, lengthmass, by = "Species")
rownames(both_phylo_test) <- both_phylo_test$Species
# Finally, repeat all of the above without the spotted porcupinefish (we will want to see
# how its exclusion as an outlier is going to change the results). This species was never
# observed on video, so there is no need to remove it from the video-only data. Note that
# the exclude() function is not safeguarded against attempts to remove taxa that are not
# present in the data frame to begin with.
nojpsp_domall <- exclude(nojp_domall, "Spotted Porcupinefish")
nojpsp_bothall <- exclude(nojp_bothall, "Spotted Porcupinefish")
nojpsp_live <- exclude(nojp_live, "Spotted Porcupinefish")
nojpsp_both <- exclude(nojp_both, "Spotted Porcupinefish")
cbi_nojpsp_domall <- cbi_hierarchy_new(nojpsp_domall)
cbi_nojpsp_bothall <- cbi_hierarchy_new(nojpsp_bothall)
cbi_nojpsp_ooolive <- cbi_hierarchy_new(nojpsp_live)
cbi_nojpsp_oooboth <- cbi_hierarchy_new(nojpsp_both)
sheet_ooo6 <- drop_taxa(sheet_ooo4, c("Cornetfish", "Juvenile Parrotfish", "Spotted Porcupinefish"))
both_ooo6 <- drop_taxa(both_ooo4, c("Cornetfish", "Juvenile Parrotfish", "Spotted Porcupinefish"))
sheet_phylo_test_nosp <- merge(cbi_nojpsp_ooolive, win_freq(sheet_ooo6), by = "Species")
sheet_phylo_test_nosp <- merge(sheet_phylo_test_nosp, lengthmass, by = "Species")
rownames(sheet_phylo_test_nosp) <- sheet_phylo_test_nosp$Species
both_phylo_test_nosp <- merge(cbi_nojpsp_oooboth, win_freq(both_ooo6), by = "Species")
both_phylo_test_nosp <- merge(both_phylo_test_nosp, lengthmass, by = "Species")
rownames(both_phylo_test_nosp) <- both_phylo_test_nosp$Species
Recalculate the properties of the one-on-one matrices created by the exclusion of juvenile parrotfish:
paste("The completeness of the live-only one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_live), "%", sep = "")
## [1] "The completeness of the live-only one-on-one matrix after the exclusion of juvenile parrotfish is 35.7142857142857%"
paste("The completeness of the video-only one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_scored), "%", sep = "")
## [1] "The completeness of the video-only one-on-one matrix after the exclusion of juvenile parrotfish is 26.6666666666667%"
paste("The completeness of the combined one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_both), "%", sep = "")
## [1] "The completeness of the combined one-on-one matrix after the exclusion of juvenile parrotfish is 34.6320346320346%"
devries(nojp_live, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.3125917
## p-value= 0.0042
devries(nojp_scored, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.3038453
## p-value= 0.07
devries(nojp_both, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2824597
## p-value= 0.0062
Now we can run phylogenetic regression:
library(nlme)
library(phytools)
## Loading required package: ape
## Loading required package: maps
fishtree <- read.tree("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/small_tree.tre")
# We will have to match the scientific names in the tip labels against the common names
# we use for the remaining data:
species_table <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", stringsAsFactors = FALSE)
fishtree$tip.label <- species_table$Common.Name[match(unlist(lapply(strsplit(fishtree$tip.label, "_"), function(x) paste(x[1], x[2], sep = " "))), species_table$Scientific.Name)]
# Make a new tree that does not include the spotted porcupinefish:
fishtree_nosp <- drop.tip(fishtree, "Spotted Porcupinefish")
# Run phylogenetic GLS on all parameter combinations (mass versus mass raised to two thirds, CBI vs. frequency of winning, spotted porcupinefish vs. no spotted porcupinefish):
pgls_CBI_Mass <- gls(CBI ~ Mass, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_CBI_Mass)
## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass
## Data: both_phylo_test
## AIC BIC logLik
## 42.96499 47.32916 -17.4825
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 1.026069
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6181935 0.3625948 1.7049155 0.1037
## Mass -0.0002343 0.0006324 -0.3705309 0.7149
##
## Correlation:
## (Intr)
## Mass -0.341
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.7446376 -0.6191146 -0.3996205 -0.2396558 3.0142764
##
## Residual standard error: 0.8129451
## Degrees of freedom: 22 total; 20 residual
pgls_CBI_Mass2_3 <- gls(CBI ~ Mass2_3, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_CBI_Mass2_3)
## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass2_3
## Data: both_phylo_test
## AIC BIC logLik
## 43.09373 47.4579 -17.54687
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 1.024889
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.5973099 0.3862164 1.5465678 0.1376
## Mass2_3 -0.0008940 0.0065095 -0.1373374 0.8921
##
## Correlation:
## (Intr)
## Mass2_3 -0.47
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.7163492 -0.5910616 -0.4011444 -0.2178379 3.0026720
##
## Residual standard error: 0.8133159
## Degrees of freedom: 22 total; 20 residual
pgls_WinFreq_Mass <- gls(WinFreq ~ Mass, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_WinFreq_Mass)
## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass
## Data: both_phylo_test
## AIC BIC logLik
## 9.610968 13.97514 -0.8054841
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7608329
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3867176 0.12406998 3.1169316 0.0054
## Mass 0.0001547 0.00024972 0.6193854 0.5427
##
## Correlation:
## (Intr)
## Mass -0.389
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.3855278 -0.9836204 -0.3401884 0.1149432 1.9053666
##
## Residual standard error: 0.2980929
## Degrees of freedom: 22 total; 20 residual
pgls_WinFreq_Mass2_3 <- gls(WinFreq ~ Mass2_3, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_WinFreq_Mass2_3)
## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass2_3
## Data: both_phylo_test
## AIC BIC logLik
## 9.341281 13.70545 -0.6706403
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7418624
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3558832 0.13389644 2.6578989 0.0151
## Mass2_3 0.0021802 0.00266542 0.8179602 0.4230
##
## Correlation:
## (Intr)
## Mass2_3 -0.551
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4406599 -0.9320170 -0.4012155 0.1434963 1.9720109
##
## Residual standard error: 0.2934742
## Degrees of freedom: 22 total; 20 residual
pgls_CBI_Mass_nosp <- gls(CBI ~ Mass, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_CBI_Mass_nosp)
## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 38.77226 42.95035 -15.38613
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.9780776
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.4131129 0.3416725 1.209090 0.2415
## Mass 0.0017121 0.0013866 1.234727 0.2320
##
## Correlation:
## (Intr)
## Mass -0.476
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.9737380 -0.6492455 -0.4673056 -0.2472223 2.9127563
##
## Residual standard error: 0.7159469
## Degrees of freedom: 21 total; 19 residual
pgls_CBI_Mass2_3_nosp <- gls(CBI ~ Mass2_3, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_CBI_Mass2_3_nosp)
## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass2_3
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 39.20199 43.38008 -15.60099
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.9938416
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.4139249 0.3738005 1.1073418 0.2820
## Mass2_3 0.0093035 0.0095417 0.9750351 0.3418
##
## Correlation:
## (Intr)
## Mass2_3 -0.551
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.8003119 -0.6502855 -0.5034040 -0.2894925 2.9422145
##
## Residual standard error: 0.7394883
## Degrees of freedom: 21 total; 19 residual
pgls_WinFreq_Mass_nosp <- gls(WinFreq ~ Mass, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_WinFreq_Mass_nosp)
## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 9.614419 13.79251 -0.8072095
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7513216
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3365943 0.13968518 2.4096639 0.0263
## Mass 0.0006368 0.00066595 0.9561862 0.3510
##
## Correlation:
## (Intr)
## Mass -0.554
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4891060 -0.9730536 -0.5088269 0.0878157 2.0020512
##
## Residual standard error: 0.2987457
## Degrees of freedom: 21 total; 19 residual
pgls_WinFreq_Mass2_3_nosp <- gls(WinFreq ~ Mass2_3, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_WinFreq_Mass2_3_nosp)
## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass2_3
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 9.600858 13.77895 -0.8004288
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7461108
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.31370935 0.15285260 2.0523652 0.0542
## Mass2_3 0.00450538 0.00465324 0.9682252 0.3451
##
## Correlation:
## (Intr)
## Mass2_3 -0.653
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.5174832 -1.0098362 -0.5279994 0.0544028 1.9944100
##
## Residual standard error: 0.2978518
## Degrees of freedom: 21 total; 19 residual
We can see that for all models involving CBI, we get very weak relationships (characterized by nearly horizontal regression lines) that are nevertheless extremely significant. A more careful examination reveals that this is because Pagel’s \(\lambda\) estimated for these models exceeds 1 (moreover, it is identical in all the four cases), even though it is only defined between 0 and 1. The cause of this is unclear (but likely related to the shape of the tree, especially to the presence of terminal branches longer than expected under the Yule process), as is the interpretability of the resulting estimate. Based on Liam Revell’s remarks in the first of the two links given below, we turn to the caper package to fit a phylogenetic GLS under the constraint that \(\lambda \leq 1\):
http://blog.phytools.org/2016/01/phylogenetic-regression-when-estimated.html http://grokbase.com/t/r/r-sig-phylo/113s5mhjgg/pagels-lambda-greater-than-one
if(!require("caper")) {install.packages("caper")}
## Loading required package: caper
## Loading required package: MASS
## Loading required package: mvtnorm
library(caper)
both_comparative <- comparative.data(data = both_phylo_test, phy = fishtree, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
both_comparative_nosp <- comparative.data(data = both_phylo_test_nosp, phy = fishtree_nosp, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = both_comparative, lambda = "ML")
summary(pgls_CBI_Mass)
##
## Call:
## pgls(formula = CBI ~ Mass, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14667 -0.02814 0.01556 0.05592 0.08199
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.010507
## upper bound : 1.000, p = 1
## 95.0% CI : (0.524, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.61306817 0.34503191 1.7768 0.09081 .
## Mass -0.00020831 0.00061615 -0.3381 0.73882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07567 on 20 degrees of freedom
## Multiple R-squared: 0.005683, Adjusted R-squared: -0.04403
## F-statistic: 0.1143 on 1 and 20 DF, p-value: 0.7388
pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = both_comparative, lambda = "ML")
summary(pgls_CBI_Mass2_3)
##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14586 -0.02786 0.01345 0.05485 0.09016
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.016049
## upper bound : 1.000, p = 1
## 95.0% CI : (0.453, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.58731317 0.37068156 1.5844 0.1288
## Mass2_3 -0.00053638 0.00644425 -0.0832 0.9345
##
## Residual standard error: 0.07587 on 20 degrees of freedom
## Multiple R-squared: 0.0003463, Adjusted R-squared: -0.04964
## F-statistic: 0.006928 on 1 and 20 DF, p-value: 0.9345
pgls_CBI_Mass_nosp <- pgls(formula = CBI ~ Mass, data = both_comparative_nosp, lambda = "ML")
summary(pgls_CBI_Mass_nosp)
##
## Call:
## pgls(formula = CBI ~ Mass, data = both_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14830 -0.04944 -0.01476 0.01939 0.18616
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.978
## lower bound : 0.000, p = 0.0081072
## upper bound : 1.000, p = 0.83889
## 95.0% CI : (0.458, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4131104 0.3416712 1.2091 0.2415
## Mass 0.0017121 0.0013866 1.2347 0.2320
##
## Residual standard error: 0.06986 on 19 degrees of freedom
## Multiple R-squared: 0.07428, Adjusted R-squared: 0.02556
## F-statistic: 1.525 on 1 and 19 DF, p-value: 0.232
pgls_CBI_Mass2_3_nosp <- pgls(formula = CBI ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
summary(pgls_CBI_Mass2_3_nosp)
##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.151213 -0.018977 0.009497 0.064289 0.139694
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.994
## lower bound : 0.000, p = 0.0090652
## upper bound : 1.000, p = 0.95565
## 95.0% CI : (0.461, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4139220 0.3737994 1.1073 0.2820
## Mass2_3 0.0093036 0.0095417 0.9750 0.3418
##
## Residual standard error: 0.07216 on 19 degrees of freedom
## Multiple R-squared: 0.04765, Adjusted R-squared: -0.00247
## F-statistic: 0.9507 on 1 and 19 DF, p-value: 0.3418
pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = both_comparative, lambda = "ML")
summary(pgls_WinFreq_Mass)
##
## Call:
## pgls(formula = WinFreq ~ Mass, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067829 -0.016253 0.002191 0.017853 0.053835
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.761
## lower bound : 0.000, p = 0.026067
## upper bound : 1.000, p = 0.060639
## 95.0% CI : (0.083, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.38671755 0.12406983 3.1169 0.005432 **
## Mass 0.00015467 0.00024972 0.6194 0.542654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02902 on 20 degrees of freedom
## Multiple R-squared: 0.01882, Adjusted R-squared: -0.03024
## F-statistic: 0.3836 on 1 and 20 DF, p-value: 0.5427
pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = both_comparative, lambda = "ML")
summary(pgls_WinFreq_Mass2_3)
##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.056556 -0.021804 -0.006782 0.014359 0.038175
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.742
## lower bound : 0.000, p = 0.028143
## upper bound : 1.000, p = 0.052267
## 95.0% CI : (0.071, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3558831 0.1338963 2.6579 0.0151 *
## Mass2_3 0.0021802 0.0026654 0.8180 0.4230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02857 on 20 degrees of freedom
## Multiple R-squared: 0.03237, Adjusted R-squared: -0.01601
## F-statistic: 0.6691 on 1 and 20 DF, p-value: 0.423
pgls_WinFreq_Mass_nosp <- pgls(formula = WinFreq ~ Mass, data = both_comparative_nosp, lambda = "ML")
summary(pgls_WinFreq_Mass_nosp)
##
## Call:
## pgls(formula = WinFreq ~ Mass, data = both_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.063430 -0.018323 -0.003256 0.012007 0.054953
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.751
## lower bound : 0.000, p = 0.022585
## upper bound : 1.000, p = 0.055727
## 95.0% CI : (0.103, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.33659421 0.13968505 2.4097 0.02627 *
## Mass 0.00063677 0.00066595 0.9562 0.35099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02915 on 19 degrees of freedom
## Multiple R-squared: 0.04591, Adjusted R-squared: -0.004304
## F-statistic: 0.9143 on 1 and 19 DF, p-value: 0.351
pgls_WinFreq_Mass2_3_nosp <- pgls(formula = WinFreq ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
summary(pgls_WinFreq_Mass2_3_nosp)
##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = both_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.061849 -0.008427 0.005194 0.019045 0.054724
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.746
## lower bound : 0.000, p = 0.023835
## upper bound : 1.000, p = 0.055536
## 95.0% CI : (0.095, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3137092 0.1528525 2.0524 0.05417 .
## Mass2_3 0.0045054 0.0046532 0.9682 0.34510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02906 on 19 degrees of freedom
## Multiple R-squared: 0.04702, Adjusted R-squared: -0.003137
## F-statistic: 0.9375 on 1 and 19 DF, p-value: 0.3451
# Check the distribution of the residuals:
par(mfrow=c(2,2))
plot(pgls_CBI_Mass)
par(mfrow=c(2,2))
plot(pgls_CBI_Mass2_3)
par(mfrow=c(2,2))
plot(pgls_CBI_Mass_nosp)
par(mfrow=c(2,2))
plot(pgls_CBI_Mass2_3_nosp)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass2_3)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass_nosp)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass2_3_nosp)
This worked great: the \(\lambda\) and p values for the winning frequency are the same as before, and while \(\lambda\) is still running up against the upper bound of 1 for CBI, the fact that it cannot exceed it results in much more sensible p values. We can now proceed to visualizing the results:
# Generate a new Figure 1:
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Figure_1.png", width = 3600, height = 3600, pointsize = 96)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.eps")
spec_dec <- function(x, k) trimws(format(round(x, k), nsmall = k)) # credit: Jeromy Anglim, Stack Overflow
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass, both_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(pgls_CBI_Mass, lwd = 5)
abline(pgls_CBI_Mass_nosp, lwd = 5, col = "gray70")
lambda <- pgls_CBI_Mass$param[2]
pseudor2 <- summary(pgls_CBI_Mass)$r.squared
pvalue <- summary(pgls_CBI_Mass)$coef[2,4]
lambda_nosp <- pgls_CBI_Mass_nosp$param[2]
pseudor2_nosp <- summary(pgls_CBI_Mass_nosp)$r.squared
pvalue_nosp <- summary(pgls_CBI_Mass_nosp)$coef[2,4]
text(600, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 1.3, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.1, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.85, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass2_3, both_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_CBI_Mass2_3, lwd = 5)
abline(pgls_CBI_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(pgls_CBI_Mass2_3)$coef[2,4]
lambda_nosp <- pgls_CBI_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(pgls_CBI_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(pgls_CBI_Mass2_3_nosp)$coef[2,4]
text(58, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 2, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 1.8, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 1.55, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass, both_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_WinFreq_Mass, lwd = 5)
abline(pgls_WinFreq_Mass_nosp, lwd = 5, col = "gray70")
lambda <- pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(pgls_WinFreq_Mass)$r.squared
pvalue <- summary(pgls_WinFreq_Mass)$coef[2,4]
lambda_nosp <- pgls_WinFreq_Mass_nosp$param[2]
pseudor2_nosp <- summary(pgls_WinFreq_Mass_nosp)$r.squared
pvalue_nosp <- summary(pgls_WinFreq_Mass_nosp)$coef[2,4]
text(600, 0.7, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 0.7 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 0.7 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 0.225, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.225 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.225 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass2_3, both_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_WinFreq_Mass2_3, lwd = 5)
abline(pgls_WinFreq_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(pgls_WinFreq_Mass2_3)$coef[2,4]
lambda_nosp <- pgls_WinFreq_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(pgls_WinFreq_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(pgls_WinFreq_Mass2_3_nosp)$coef[2,4]
text(58, 0.98, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 0.98 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 0.98 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.225, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.225 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.225 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()
Below, we conduct another test in which the \(\alpha\) coefficients (individual fighting abilities) are all set to 1, and therefore effectively disregarded. We also evaluate the fit of another model, which includes body size as the only parameter. Note that neither model requires us to repeat the test for body mass raised to 2/3. The group size-only model does not include this parameter at all, while the body size-only model will always yield the same predictions under both mass scalings, since \(m^{\frac{2}{3}}\) is a monotonically increasing function of \(m\).
# Add a new column to both_group with predictions based on group size alone:
for(i in 1:nrow(both_group)) {
if (both_group$Group.A[i] > both_group$Group.B[i]) {
both_group[i,14] <- both_group$Species.A[i]
} else {
both_group[i,14] <- both_group$Species.B[i]
}
}
# Add a new column to both_group with predictions based on body size alone:
for(i in 1:nrow(both_group)) {
if (both_group$BodySize.A[i] > both_group$BodySize.B[i]) {
both_group[i,15] <- both_group$Species.A[i]
} else {
both_group[i,15] <- both_group$Species.B[i]
}
}
colnames(both_group)[14:15] <- c("Group.Size.Only.Prediction","Body.Size.Only.Prediction")
# Test if the group size-only model works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Group.Size.Only.Prediction[i])) == TRUE) {
both_group[i,16] <- 1
} else {
both_group[i,16] <- 0
}
}
# Test if the body size-only model works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Body.Size.Only.Prediction[i])) == TRUE) {
both_group[i,17] <- 1
} else {
both_group[i,17] <- 0
}
}
colnames(both_group)[16:17] <- c("Group.Size.Only.True", "Body.Size.Only.True")
# Proportion of successful predictions:
sum(both_group$Group.Size.Only.True)/nrow(both_group)
## [1] 0.2289157
sum(both_group$Body.Size.Only.True)/nrow(both_group)
## [1] 0.6987952
# Test significance:
dbinom(sum(both_group$Group.Size.Only.True), size = nrow(both_group), prob = 1/2)
## [1] 2.643039e-07
dbinom(sum(both_group$Body.Size.Only.True), size = nrow(both_group), prob = 1/2)
## [1] 0.0001118917
In what cases did the body size-only model work but the next best model (Lanchester’s linear law) did not?
both_group[both_group$Body.Size.Only.True == 1 & both_group$Linear.True == 0,]
## Species.A Species.B Group.A Group.B
## 59 Damselfish Scissortail Sergeant 1 2
## 87 Damselfish Sixbar Wrasse 1 2
## 246 Damselfish Juvenile Parrotfish 1 20
## 259 Regal Angelfish Brown Tang 1 2
## 273 Orange-lined Triggerfish Threeband Pennantfish 1 2
## 276 Orange-lined Triggerfish Brown Tang 1 2
## 324 Damselfish Scissortail Sergeant 1 3
## 329 Damselfish Sixbar Wrasse 1 2
## 487 Damselfish Juvenile Parrotfish 1 7
## 515 Damselfish Juvenile Parrotfish 1 14
## 554 Damselfish Scissortail Sergeant 1 2
## Winner BodySize.A BodySize.B Linear.Law.Prediction
## 59 Damselfish 38.88513 25.35153 Scissortail Sergeant
## 87 Damselfish 38.88513 22.72689 Sixbar Wrasse
## 246 Damselfish 38.88513 7.49326 Juvenile Parrotfish
## 259 Regal Angelfish 128.17052 78.11946 Brown Tang
## 273 Orange-lined Triggerfish 108.96596 87.65761 Threeband Pennantfish
## 276 Orange-lined Triggerfish 108.96596 78.11946 Brown Tang
## 324 Damselfish 38.88513 25.35153 Scissortail Sergeant
## 329 Damselfish 38.88513 22.72689 Sixbar Wrasse
## 487 Damselfish 38.88513 7.49326 Juvenile Parrotfish
## 515 Damselfish 38.88513 7.49326 Juvenile Parrotfish
## 554 Damselfish 38.88513 25.35153 Scissortail Sergeant
## Square.Law.Prediction Linear.True Square.True
## 59 Scissortail Sergeant 0 0
## 87 Sixbar Wrasse 0 0
## 246 Juvenile Parrotfish 0 0
## 259 Brown Tang 0 0
## 273 Threeband Pennantfish 0 0
## 276 Brown Tang 0 0
## 324 Scissortail Sergeant 0 0
## 329 Sixbar Wrasse 0 0
## 487 Juvenile Parrotfish 0 0
## 515 Juvenile Parrotfish 0 0
## 554 Scissortail Sergeant 0 0
## Diet.A Diet.B
## 59 Omnivore Omnivore
## 87 Omnivore Benthic Invertebrate Consumer
## 246 Omnivore Herbivore/Detritivore
## 259 Benthic Invertebrate Consumer Herbivore/Detritivore
## 273 Omnivore Corallivore
## 276 Omnivore Herbivore/Detritivore
## 324 Omnivore Omnivore
## 329 Omnivore Benthic Invertebrate Consumer
## 487 Omnivore Herbivore/Detritivore
## 515 Omnivore Herbivore/Detritivore
## 554 Omnivore Omnivore
## Group.Size.Only.Prediction Body.Size.Only.Prediction
## 59 Scissortail Sergeant Damselfish
## 87 Sixbar Wrasse Damselfish
## 246 Juvenile Parrotfish Damselfish
## 259 Brown Tang Regal Angelfish
## 273 Threeband Pennantfish Orange-lined Triggerfish
## 276 Brown Tang Orange-lined Triggerfish
## 324 Scissortail Sergeant Damselfish
## 329 Sixbar Wrasse Damselfish
## 487 Juvenile Parrotfish Damselfish
## 515 Juvenile Parrotfish Damselfish
## 554 Scissortail Sergeant Damselfish
## Group.Size.Only.True Body.Size.Only.True
## 59 0 1
## 87 0 1
## 246 0 1
## 259 0 1
## 273 0 1
## 276 0 1
## 324 0 1
## 329 0 1
## 487 0 1
## 515 0 1
## 554 0 1
8 out of the 11 interactions in this category involved damselfish. What happens if we remove them from the dataset?
both_group3 <- both_group[both_group$Species.A != "Damselfish" & both_group$Species.B != "Damselfish", ]
nrow(both_group3)
## [1] 43
# Proportion of successful predictions:
sum(both_group3$Linear.True)/nrow(both_group3)
## [1] 0.8372093
sum(both_group3$Square.True)/nrow(both_group3)
## [1] 0.7209302
sum(both_group3$Group.Size.Only.True)/nrow(both_group3)
## [1] 0.255814
sum(both_group3$Body.Size.Only.True)/nrow(both_group3)
## [1] 0.7906977
# Test significance:
dbinom(sum(both_group3$Linear.True), size = nrow(both_group3), prob = 1/2)
## [1] 3.663458e-06
dbinom(sum(both_group3$Square.True), size = nrow(both_group3), prob = 1/2)
## [1] 0.001743806
dbinom(sum(both_group3$Group.Size.Only.True), size = nrow(both_group3), prob = 1/2)
## [1] 0.0006539272
dbinom(sum(both_group3$Body.Size.Only.True), size = nrow(both_group3), prob = 1/2)
## [1] 6.411051e-05
# With mass raised to the two-thirds power. Remember, this can change the results for the
# linear and square laws, but not for the body size-only or group size-only models.
both_group4 <- both_group2[both_group2$Species.A != "Damselfish" & both_group2$Species.B != "Damselfish" ,]
nrow(both_group4)
## [1] 43
# Proportion of successful predictions:
sum(both_group4$Linear.True)/nrow(both_group4)
## [1] 0.7906977
sum(both_group4$Square.True)/nrow(both_group4)
## [1] 0.5348837
# Test significance:
dbinom(sum(both_group4$Linear.True), size = nrow(both_group4), prob = 1/2)
## [1] 6.411051e-05
dbinom(sum(both_group4$Square.True), size = nrow(both_group4), prob = 1/2)
## [1] 0.1092038
In the next step, we calculate the statistical power of the linear law, square law, and body size-only model tests. All of the below is based on http://www.calvin.edu/~scofield/courses/m343/F15/handouts/binomialTestPower.pdf
# First, find the rejection region for n = 83 (the number of data points in both_group)
# under the null assumption (probability of successful prediction = 0.5) and significance
# level alpha = 0.05:
qbinom(0.025, nrow(both_group), 0.5) # 33
## [1] 33
qbinom(0.975, nrow(both_group), 0.5) # 50
## [1] 50
# Since P(33 \leq x \geq 50) is probably not equal to exactly 0.95, we need to find out
# if the sum of the values below is larger or smaller than this significance threshold:
pbinom(33, nrow(both_group), 0.5) + (1 - pbinom(50, nrow(both_group), 0.5))
## [1] 0.06297226
pbinom(32, nrow(both_group), 0.5) + (1 - pbinom(51, nrow(both_group), 0.5))
## [1] 0.03752953
# Let's use 32 and 51 as the boundaries of the rejection region.
# Compute power as 1 - beta:
1 - (pbinom(51, nrow(both_group), sum(both_group$Body.Size.Only.True)/nrow(both_group)) -
pbinom(32, nrow(both_group), sum(both_group$Body.Size.Only.True)/nrow(both_group)))
## [1] 0.9377759
Compute the maximum likelihood estimate under the binomial and use it to compute AIC scores:
# -2*log(L^hat) for the body size-only model:
-2*log(dbinom(sum(both_group$Body.Size.Only.True),
nrow(both_group),
sum(both_group$Body.Size.Only.True)/nrow(both_group)))
## [1] 4.705887
# -2*log(L^hat) for the linear law:
-2*log(dbinom(sum(both_group$Linear.True),
nrow(both_group),
sum(both_group$Linear.True)/nrow(both_group)))
## [1] 4.765549
Finally, see how the body size-only and group-size only models perform when fitted to the datasets subsampled by diet:
test_omni
## Species.A Species.B Group.A Group.B
## 59 Damselfish Scissortail Sergeant 1 2
## 113 Moorish Idol Lemonpeel Angelfish 1 2
## 132 Damselfish Lemonpeel Angelfish 1 2
## 187 Dusky Angelfish Scissortail Sergeant 1 3
## 288 Orange-lined Triggerfish Lemonpeel Angelfish 1 2
## 324 Damselfish Scissortail Sergeant 1 3
## 353 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 355 Damselfish Lemonpeel Angelfish 1 3
## 357 Damselfish Lemonpeel Angelfish 1 3
## 358 Damselfish Lemonpeel Angelfish 2 1
## 363 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 366 Moorish Idol Lemonpeel Angelfish 1 3
## 368 Damselfish Lemonpeel Angelfish 1 3
## 373 Damselfish Lemonpeel Angelfish 1 3
## 440 Damselfish Lemonpeel Angelfish 1 2
## 554 Damselfish Scissortail Sergeant 1 2
## Winner BodySize.A BodySize.B
## 59 Damselfish 38.88513 25.351526
## 113 Moorish Idol 171.38610 7.490788
## 132 Damselfish 38.88513 7.490788
## 187 Scissortail Sergeant 10.99104 25.351526
## 288 Orange-lined Triggerfish 108.96596 7.490788
## 324 Damselfish 38.88513 25.351526
## 353 Orange-lined Triggerfish 108.96596 7.490788
## 355 Damselfish 38.88513 7.490788
## 357 Damselfish 38.88513 7.490788
## 358 Damselfish 38.88513 7.490788
## 363 Orange-lined Triggerfish 108.96596 7.490788
## 366 Moorish Idol 171.38610 7.490788
## 368 Damselfish 38.88513 7.490788
## 373 Damselfish 38.88513 7.490788
## 440 Damselfish 38.88513 7.490788
## 554 Damselfish 38.88513 25.351526
## Linear.Law.Prediction Square.Law.Prediction Linear.True
## 59 Scissortail Sergeant Scissortail Sergeant 0
## 113 Moorish Idol Moorish Idol 1
## 132 Damselfish Damselfish 1
## 187 Scissortail Sergeant Scissortail Sergeant 1
## 288 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 324 Scissortail Sergeant Scissortail Sergeant 0
## 353 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 355 Damselfish Lemonpeel Angelfish 1
## 357 Damselfish Lemonpeel Angelfish 1
## 358 Damselfish Damselfish 1
## 363 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 366 Moorish Idol Moorish Idol 1
## 368 Damselfish Lemonpeel Angelfish 1
## 373 Damselfish Lemonpeel Angelfish 1
## 440 Damselfish Damselfish 1
## 554 Scissortail Sergeant Scissortail Sergeant 0
## Square.True Diet.A Diet.B
## 59 0 Omnivore Omnivore
## 113 1 Omnivore Omnivore
## 132 1 Omnivore Omnivore
## 187 1 Omnivore Omnivore
## 288 1 Omnivore Omnivore
## 324 0 Omnivore Omnivore
## 353 1 Omnivore Omnivore
## 355 0 Omnivore Omnivore
## 357 0 Omnivore Omnivore
## 358 1 Omnivore Omnivore
## 363 1 Omnivore Omnivore
## 366 1 Omnivore Omnivore
## 368 0 Omnivore Omnivore
## 373 0 Omnivore Omnivore
## 440 1 Omnivore Omnivore
## 554 0 Omnivore Omnivore
test_omni2
## Species.A Species.B Group.A Group.B
## 59 Damselfish Scissortail Sergeant 1 2
## 113 Moorish Idol Lemonpeel Angelfish 1 2
## 132 Damselfish Lemonpeel Angelfish 1 2
## 187 Dusky Angelfish Scissortail Sergeant 1 3
## 288 Orange-lined Triggerfish Lemonpeel Angelfish 1 2
## 324 Damselfish Scissortail Sergeant 1 3
## 353 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 355 Damselfish Lemonpeel Angelfish 1 3
## 357 Damselfish Lemonpeel Angelfish 1 3
## 358 Damselfish Lemonpeel Angelfish 2 1
## 363 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 366 Moorish Idol Lemonpeel Angelfish 1 3
## 368 Damselfish Lemonpeel Angelfish 1 3
## 373 Damselfish Lemonpeel Angelfish 1 3
## 440 Damselfish Lemonpeel Angelfish 1 2
## 554 Damselfish Scissortail Sergeant 1 2
## Winner BodySize.A BodySize.B
## 59 Damselfish 11.477723 8.629840
## 113 Moorish Idol 30.854397 3.828409
## 132 Damselfish 11.477723 3.828409
## 187 Scissortail Sergeant 4.943401 8.629840
## 288 Orange-lined Triggerfish 22.813604 3.828409
## 324 Damselfish 11.477723 8.629840
## 353 Orange-lined Triggerfish 22.813604 3.828409
## 355 Damselfish 11.477723 3.828409
## 357 Damselfish 11.477723 3.828409
## 358 Damselfish 11.477723 3.828409
## 363 Orange-lined Triggerfish 22.813604 3.828409
## 366 Moorish Idol 30.854397 3.828409
## 368 Damselfish 11.477723 3.828409
## 373 Damselfish 11.477723 3.828409
## 440 Damselfish 11.477723 3.828409
## 554 Damselfish 11.477723 8.629840
## Linear.Law.Prediction Square.Law.Prediction Linear.True
## 59 Scissortail Sergeant Scissortail Sergeant 0
## 113 Moorish Idol Moorish Idol 1
## 132 Damselfish Lemonpeel Angelfish 1
## 187 Scissortail Sergeant Scissortail Sergeant 1
## 288 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 324 Scissortail Sergeant Scissortail Sergeant 0
## 353 Orange-lined Triggerfish Lemonpeel Angelfish 1
## 355 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 357 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 358 Damselfish Damselfish 1
## 363 Orange-lined Triggerfish Lemonpeel Angelfish 1
## 366 Moorish Idol Lemonpeel Angelfish 1
## 368 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 373 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 440 Damselfish Lemonpeel Angelfish 1
## 554 Scissortail Sergeant Scissortail Sergeant 0
## Square.True Diet.A Diet.B
## 59 0 Omnivore Omnivore
## 113 1 Omnivore Omnivore
## 132 0 Omnivore Omnivore
## 187 1 Omnivore Omnivore
## 288 1 Omnivore Omnivore
## 324 0 Omnivore Omnivore
## 353 0 Omnivore Omnivore
## 355 0 Omnivore Omnivore
## 357 0 Omnivore Omnivore
## 358 1 Omnivore Omnivore
## 363 0 Omnivore Omnivore
## 366 0 Omnivore Omnivore
## 368 0 Omnivore Omnivore
## 373 0 Omnivore Omnivore
## 440 0 Omnivore Omnivore
## 554 0 Omnivore Omnivore
test_herb
## Species.A Species.B Group.A Group.B
## 116 Brown Tang Juvenile Parrotfish 1 4
## 182 Dark-capped Parrotfish Juvenile Parrotfish 1 10
## 264 Bullethead Parrotfish Brown Tang 1 2
## 282 Dark-capped Parrotfish Juvenile Parrotfish 1 20
## Winner BodySize.A BodySize.B Linear.Law.Prediction
## 116 Brown Tang 78.11946 7.49326 Brown Tang
## 182 Dark-capped Parrotfish 311.91661 7.49326 Dark-capped Parrotfish
## 264 Bullethead Parrotfish 339.38694 78.11946 Bullethead Parrotfish
## 282 Dark-capped Parrotfish 311.91661 7.49326 Dark-capped Parrotfish
## Square.Law.Prediction Linear.True Square.True Diet.A
## 116 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 182 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 264 Bullethead Parrotfish 1 1 Herbivore/Detritivore
## 282 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## Diet.B
## 116 Herbivore/Detritivore
## 182 Herbivore/Detritivore
## 264 Herbivore/Detritivore
## 282 Herbivore/Detritivore
test_herb2
## Species.A Species.B Group.A Group.B
## 116 Brown Tang Juvenile Parrotfish 1 4
## 182 Dark-capped Parrotfish Juvenile Parrotfish 1 10
## 264 Bullethead Parrotfish Brown Tang 1 2
## 282 Dark-capped Parrotfish Juvenile Parrotfish 1 20
## Winner BodySize.A BodySize.B Linear.Law.Prediction
## 116 Brown Tang 18.27425 3.829251 Brown Tang
## 182 Dark-capped Parrotfish 45.99306 3.829251 Dark-capped Parrotfish
## 264 Bullethead Parrotfish 48.65529 18.274246 Bullethead Parrotfish
## 282 Dark-capped Parrotfish 45.99306 3.829251 Juvenile Parrotfish
## Square.Law.Prediction Linear.True Square.True Diet.A
## 116 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 182 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 264 Brown Tang 1 0 Herbivore/Detritivore
## 282 Juvenile Parrotfish 0 0 Herbivore/Detritivore
## Diet.B
## 116 Herbivore/Detritivore
## 182 Herbivore/Detritivore
## 264 Herbivore/Detritivore
## 282 Herbivore/Detritivore
# The function below adds columns "Group.Size.Only.Prediction",
# "Body.Size.Only.Prediction", "Group.Size.Only.True", and "Body.Size.Only.True" to an
# arbitrary data frame containing columns named "Group.A", "Group.B", "BodySizeA",
# "BodySizeB", and "Winner":
test_extra_two_models <- function(dataframe) {
n <- ncol(dataframe)
for(i in 1:nrow(dataframe)) {
if (dataframe$Group.A[i] > dataframe$Group.B[i]) {
dataframe[i, n+1] <- dataframe$Species.A[i]
} else {
dataframe[i, n+1] <- dataframe$Species.B[i]
}
if (dataframe$BodySize.A[i] > dataframe$BodySize.B[i]) {
dataframe[i, n+2] <- dataframe$Species.A[i]
} else {
dataframe[i, n+2] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+1])) == TRUE) {
dataframe[i, n+3] <- 1
} else {
dataframe[i, n+3] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+2])) == TRUE) {
dataframe[i, n+4] <- 1
} else {
dataframe[i, n+4] <- 0
}
}
colnames(dataframe)[(n+1):(n+4)] <- c("Group.Size.Only.Prediction",
"Body.Size.Only.Prediction",
"Group.Size.Only.True",
"Body.Size.Only.True")
return(dataframe)
}
test_omni <- test_extra_two_models(test_omni)
test_herb <- test_extra_two_models(test_herb)
# Proportion of successful predictions:
sum(test_omni$Group.Size.Only.True)/nrow(test_omni)
## [1] 0.125
sum(test_omni$Body.Size.Only.True)/nrow(test_omni)
## [1] 1
sum(test_herb$Group.Size.Only.True)/nrow(test_herb)
## [1] 0
sum(test_herb$Body.Size.Only.True)/nrow(test_herb)
## [1] 1
# Test significance:
dbinom(sum(test_omni$Group.Size.Only.True), size = nrow(test_omni), prob = 1/2)
## [1] 0.001831055
dbinom(sum(test_omni$Body.Size.Only.True), size = nrow(test_omni), prob = 1/2)
## [1] 1.525879e-05
dbinom(sum(test_herb$Group.Size.Only.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.0625
dbinom(sum(test_herb$Body.Size.Only.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.0625
We have already alculated the properties of matrices based on live-only and video-only observations. In the next step, we compare the hierarchies:
# Subsample each taxon-by-taxon dominance matrix down to the shared set of 15 species
nojp_live_shared <- exclude(nojp_live, setdiff(rownames(nojp_live), rownames(nojp_scored)))
nojp_scored_shared <- exclude(nojp_scored, setdiff(rownames(nojp_scored), rownames(nojp_live)))
nojp_both_shared <- exclude(nojp_both, rownames(nojp_both)[!rownames(nojp_both) %in% intersect(rownames(nojp_live), rownames(nojp_scored))])
# Recalculate the CBI for the 15 species based on the live-only and video-only datasets:
live_subsampled_cbi <- cbi_hierarchy_new(nojp_live_shared)
colnames(live_subsampled_cbi)[2] <- "Live"
video_subsampled_cbi <- cbi_hierarchy_new(nojp_scored_shared)
colnames(video_subsampled_cbi)[2] <- "Video"
subsampled_cbi <- merge(live_subsampled_cbi, video_subsampled_cbi, by = "Species")
both_subsampled_cbi <- cbi_hierarchy_new(nojp_both_shared)
subsampled_cbi <- merge(subsampled_cbi, both_subsampled_cbi, by = "Species")
colnames(subsampled_cbi)[4] <- "Both"
# Recalculate the frequency of winning for the 15 species based on the live-only and
# video-only datasets:
live_subsampled_winfreq <- win_freq(sheet_ooo5)[win_freq(sheet_ooo5)$Species %in% intersect(win_freq(sheet_ooo5)$Species, win_freq(scored_ooo5)$Species),]
live_subsampled_winfreq <- live_subsampled_winfreq[order(live_subsampled_winfreq$Species),]
colnames(live_subsampled_winfreq)[2] <- "Live"
video_subsampled_winfreq <- win_freq(scored_ooo5)[win_freq(scored_ooo5)$Species %in% intersect(win_freq(sheet_ooo5)$Species, win_freq(scored_ooo5)$Species),]
video_subsampled_winfreq <- video_subsampled_winfreq[order(video_subsampled_winfreq$Species),]
colnames(video_subsampled_winfreq)[2] <- "Video"
subsampled_winfreq <- merge(live_subsampled_winfreq, video_subsampled_winfreq, by = "Species")
subsampled_winfreq <- merge(subsampled_winfreq, win_freq(both_ooo5), by = "Species")
colnames(subsampled_winfreq)[4] <- "Both"
if(!require("RcmdrMisc")) {install.packages("RcmdrMisc")} # Beware: has many dependencies!
## Loading required package: RcmdrMisc
## Loading required package: car
## Loading required package: carData
## Loading required package: sandwich
library(RcmdrMisc)
rcorr.adjust(as.matrix(subsampled_cbi[,-1]), type = "spearman")
##
## Spearman correlations:
## Live Video Both
## Live 1.0000 0.2898 0.9357
## Video 0.2898 1.0000 0.3792
## Both 0.9357 0.3792 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.2948 <.0001
## Video 0.2948 0.1633
## Both <.0001 0.1633
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.3266 <.0001
## Video 0.3266 0.3266
## Both <.0001 0.3266
rcorr.adjust(as.matrix(subsampled_cbi[,-1]), type = "pearson")
##
## Pearson correlations:
## Live Video Both
## Live 1.0000 0.6703 0.8643
## Video 0.6703 1.0000 0.8914
## Both 0.8643 0.8914 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.0062 <.0001
## Video 0.0062 <.0001
## Both <.0001 <.0001
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.0062 <.0001
## Video 0.0062 <.0001
## Both <.0001 <.0001
rcorr.adjust(as.matrix(subsampled_winfreq[,-1]), type = "spearman")
##
## Spearman correlations:
## Live Video Both
## Live 1.0000 0.6188 0.8959
## Video 0.6188 1.0000 0.8307
## Both 0.8959 0.8307 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.0139 <.0001
## Video 0.0139 0.0001
## Both <.0001 0.0001
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.0139 <.0001
## Video 0.0139 0.0003
## Both <.0001 0.0003
rcorr.adjust(as.matrix(subsampled_winfreq[,-1]), type = "pearson")
##
## Pearson correlations:
## Live Video Both
## Live 1.0000 0.6208 0.9133
## Video 0.6208 1.0000 0.8731
## Both 0.9133 0.8731 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.0135 <.0001
## Video 0.0135 <.0001
## Both <.0001 <.0001
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.0135 <.0001
## Video 0.0135 <.0001
## Both <.0001 <.0001
Re-test the relationship between body size and individual fighting ability using live-only and video-only data:
# Appropriately subsample 'fishtree':
fishtree_live <- drop.tip(fishtree, fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_live)])
fishtree_live_nosp <- drop.tip(fishtree, c("Spotted Porcupinefish", fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_live)]))
fishtree_video <- drop.tip(fishtree, fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_scored)])
library(caper)
live_comparative <- comparative.data(data = sheet_phylo_test, phy = fishtree_live, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
live_comparative_nosp <- comparative.data(data = sheet_phylo_test_nosp, phy = fishtree_live_nosp, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
video_comparative <- comparative.data(data = scored_phylo_test, phy = fishtree_video, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
live_pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = live_comparative, lambda = "ML")
summary(live_pgls_CBI_Mass)
##
## Call:
## pgls(formula = CBI ~ Mass, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12134 -0.09154 -0.02659 0.07927 0.14682
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00051509
## upper bound : 1.000, p = 1
## 95.0% CI : (0.803, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.04492043 0.41860472 2.4962 0.02192 *
## Mass -0.00048654 0.00074678 -0.6515 0.52251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09167 on 19 degrees of freedom
## Multiple R-squared: 0.02185, Adjusted R-squared: -0.02963
## F-statistic: 0.4245 on 1 and 19 DF, p-value: 0.5225
live_pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = live_comparative, lambda = "ML"); summary(live_pgls_CBI_Mass2_3)
##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14217 -0.08418 -0.02553 0.07819 0.14423
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00083993
## upper bound : 1.000, p = 1
## 95.0% CI : (0.792, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0018370 0.4523589 2.2147 0.0392 *
## Mass2_3 -0.0018632 0.0078794 -0.2365 0.8156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09256 on 19 degrees of freedom
## Multiple R-squared: 0.002934, Adjusted R-squared: -0.04954
## F-statistic: 0.05591 on 1 and 19 DF, p-value: 0.8156
live_pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = live_comparative, lambda = "ML"); summary(live_pgls_WinFreq_Mass)
##
## Call:
## pgls(formula = WinFreq ~ Mass, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.037100 -0.016062 0.006017 0.029263 0.070629
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.663
## lower bound : 0.000, p = 0.10993
## upper bound : 1.000, p = 0.01489
## 95.0% CI : (NA, 0.973)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.41376043 0.13109752 3.1561 0.005201 **
## Mass 0.00020690 0.00027338 0.7568 0.458435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03139 on 19 degrees of freedom
## Multiple R-squared: 0.02926, Adjusted R-squared: -0.02183
## F-statistic: 0.5728 on 1 and 19 DF, p-value: 0.4584
live_pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = live_comparative, lambda = "ML"); summary(live_pgls_WinFreq_Mass2_3)
##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067679 -0.011879 -0.003656 0.017396 0.048635
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.631
## lower bound : 0.000, p = 0.13084
## upper bound : 1.000, p = 0.011061
## 95.0% CI : (NA, 0.966)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3593571 0.1397527 2.5714 0.01869 *
## Mass2_3 0.0033572 0.0028881 1.1624 0.25946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03047 on 19 degrees of freedom
## Multiple R-squared: 0.06639, Adjusted R-squared: 0.01726
## F-statistic: 1.351 on 1 and 19 DF, p-value: 0.2595
live_pgls_CBI_Mass_nosp <- pgls(formula = CBI ~ Mass, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_CBI_Mass_nosp)
##
## Call:
## pgls(formula = CBI ~ Mass, data = live_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140248 -0.004507 0.026080 0.089520 0.192394
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00037668
## upper bound : 1.000, p = 1
## 95.0% CI : (0.787, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8128104 0.4256586 1.9095 0.07226 .
## Mass 0.0017328 0.0016737 1.0353 0.31424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0873 on 18 degrees of freedom
## Multiple R-squared: 0.0562, Adjusted R-squared: 0.003767
## F-statistic: 1.072 on 1 and 18 DF, p-value: 0.3142
live_pgls_CBI_Mass2_3_nosp <- pgls(formula = CBI ~ Mass2_3, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_CBI_Mass2_3_nosp)
##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = live_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140775 -0.004609 0.022412 0.089080 0.188876
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00039991
## upper bound : 1.000, p = 1
## 95.0% CI : (0.785, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.750749 0.450167 1.6677 0.1127
## Mass2_3 0.012323 0.011456 1.0756 0.2963
##
## Residual standard error: 0.0871 on 18 degrees of freedom
## Multiple R-squared: 0.06039, Adjusted R-squared: 0.008192
## F-statistic: 1.157 on 1 and 18 DF, p-value: 0.2963
live_pgls_WinFreq_Mass_nosp <- pgls(formula = WinFreq ~ Mass, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_WinFreq_Mass_nosp)
##
## Call:
## pgls(formula = WinFreq ~ Mass, data = live_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.014123 0.002685 0.014639 0.036257 0.064441
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.647
## lower bound : 0.000, p = 0.079667
## upper bound : 1.000, p = 0.0081626
## 95.0% CI : (NA, 0.959)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.31444575 0.14286027 2.2011 0.04102 *
## Mass 0.00116556 0.00070943 1.6430 0.11775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03031 on 18 degrees of freedom
## Multiple R-squared: 0.1304, Adjusted R-squared: 0.0821
## F-statistic: 2.699 on 1 and 18 DF, p-value: 0.1177
live_pgls_WinFreq_Mass2_3_nosp <- pgls(formula = WinFreq ~ Mass2_3, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_WinFreq_Mass2_3_nosp)
##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = live_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.061269 -0.013488 -0.002605 0.015575 0.042127
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.643
## lower bound : 0.000, p = 0.079968
## upper bound : 1.000, p = 0.0076034
## 95.0% CI : (NA, 0.957)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2629693 0.1564604 1.6807 0.11008
## Mass2_3 0.0087316 0.0049639 1.7590 0.09556 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02997 on 18 degrees of freedom
## Multiple R-squared: 0.1467, Adjusted R-squared: 0.09928
## F-statistic: 3.094 on 1 and 18 DF, p-value: 0.09556
video_pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = video_comparative, lambda = "ML")
summary(video_pgls_CBI_Mass)
##
## Call:
## pgls(formula = CBI ~ Mass, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12487 -0.07136 -0.04121 0.01673 0.10692
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.049
## lower bound : 0.000, p = 0.87521
## upper bound : 1.000, p = 0.00084113
## 95.0% CI : (NA, 0.719)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6174704 0.2932005 2.1060 0.05374 .
## Mass -0.0014836 0.0020582 -0.7208 0.48287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07338 on 14 degrees of freedom
## Multiple R-squared: 0.03579, Adjusted R-squared: -0.03308
## F-statistic: 0.5196 on 1 and 14 DF, p-value: 0.4829
video_pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = video_comparative, lambda = "ML"); summary(video_pgls_CBI_Mass2_3)
##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11810 -0.05457 -0.01404 0.03946 0.10942
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.054
## lower bound : 0.000, p = 0.86644
## upper bound : 1.000, p = 0.00091886
## 95.0% CI : (NA, 0.725)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6628957 0.3641409 1.8204 0.09013 .
## Mass2_3 -0.0099129 0.0152288 -0.6509 0.52563
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07362 on 14 degrees of freedom
## Multiple R-squared: 0.02938, Adjusted R-squared: -0.03995
## F-statistic: 0.4237 on 1 and 14 DF, p-value: 0.5256
video_pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = video_comparative, lambda = "ML"); summary(video_pgls_WinFreq_Mass)
##
## Call:
## pgls(formula = WinFreq ~ Mass, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.035535 -0.011008 0.001371 0.023994 0.063956
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.536
## lower bound : 0.000, p = 0.21789
## upper bound : 1.000, p = 0.15603
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28375197 0.13948795 2.0342 0.06133 .
## Mass -0.00043991 0.00083487 -0.5269 0.60650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02897 on 14 degrees of freedom
## Multiple R-squared: 0.01945, Adjusted R-squared: -0.05059
## F-statistic: 0.2776 on 1 and 14 DF, p-value: 0.6065
video_pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = video_comparative, lambda = "ML"); summary(video_pgls_WinFreq_Mass2_3)
##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.031340 -0.010905 -0.002889 0.022627 0.064360
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.558
## lower bound : 0.000, p = 0.21225
## upper bound : 1.000, p = 0.16822
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2890066 0.1631528 1.7714 0.09825 .
## Mass2_3 -0.0025029 0.0059583 -0.4201 0.68081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02923 on 14 degrees of freedom
## Multiple R-squared: 0.01245, Adjusted R-squared: -0.05809
## F-statistic: 0.1765 on 1 and 14 DF, p-value: 0.6808
# Check the distribution of the residuals:
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass)
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass2_3)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass2_3)
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass_nosp)
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass2_3_nosp)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass_nosp)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass2_3_nosp)
par(mfrow=c(2,2))
plot(video_pgls_CBI_Mass)
par(mfrow=c(2,2))
plot(video_pgls_CBI_Mass2_3)
par(mfrow=c(2,2))
plot(video_pgls_WinFreq_Mass)
par(mfrow=c(2,2))
plot(video_pgls_WinFreq_Mass2_3)
Generate dataset-specific versions of Figure 1:
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Supp_Fig1_live.png", width = 3600, height = 3600, pointsize = 96)
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass, sheet_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(live_pgls_CBI_Mass, lwd = 5)
abline(live_pgls_CBI_Mass_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_CBI_Mass$param[2]
pseudor2 <- summary(live_pgls_CBI_Mass)$r.squared
pvalue <- summary(live_pgls_CBI_Mass)$coef[2,4]
lambda_nosp <- live_pgls_CBI_Mass_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_CBI_Mass_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_CBI_Mass_nosp)$coef[2,4]
text(600, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 1.6, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.15, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass2_3, sheet_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_CBI_Mass2_3, lwd = 5)
abline(live_pgls_CBI_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(live_pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(live_pgls_CBI_Mass2_3)$coef[2,4]
lambda_nosp <- live_pgls_CBI_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_CBI_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_CBI_Mass2_3_nosp)$coef[2,4]
text(58, 2.75, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 2.55, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 2.3, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.6, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.15, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass, sheet_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_WinFreq_Mass, lwd = 5)
abline(live_pgls_WinFreq_Mass_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(live_pgls_WinFreq_Mass)$r.squared
pvalue <- summary(live_pgls_WinFreq_Mass)$coef[2,4]
lambda_nosp <- live_pgls_WinFreq_Mass_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_WinFreq_Mass_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_WinFreq_Mass_nosp)$coef[2,4]
text(600, 0.9, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 0.9 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 0.9 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 0.35, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.35 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.35 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass2_3, sheet_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_WinFreq_Mass2_3, lwd = 5)
abline(live_pgls_WinFreq_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(live_pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(live_pgls_WinFreq_Mass2_3)$coef[2,4]
lambda_nosp <- live_pgls_WinFreq_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_WinFreq_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_WinFreq_Mass2_3_nosp)$coef[2,4]
text(58, 0.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 0.5 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 0.5 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.2, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.2 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.2 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Supp_Fig1_video.png", width = 3600, height = 3600, pointsize = 96)
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass, scored_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = "black")
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(video_pgls_CBI_Mass, lwd = 5)
lambda <- video_pgls_CBI_Mass$param[2]
pseudor2 <- summary(video_pgls_CBI_Mass)$r.squared
pvalue <- summary(video_pgls_CBI_Mass)$coef[2,4]
text(150, 2.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(150, 2.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(150, 2.05, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass2_3, scored_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = "black")
abline(video_pgls_CBI_Mass2_3, lwd = 5)
lambda <- video_pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(video_pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(video_pgls_CBI_Mass2_3)$coef[2,4]
text(25, 2.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(25, 2.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(25, 2.05, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass, scored_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = "black")
abline(video_pgls_WinFreq_Mass, lwd = 5)
lambda <- video_pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(video_pgls_WinFreq_Mass)$r.squared
pvalue <- summary(video_pgls_WinFreq_Mass)$coef[2,4]
text(150, 0.55, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(150, 0.55 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(150, 0.55 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass2_3, scored_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = "black")
abline(video_pgls_WinFreq_Mass2_3, lwd = 5)
lambda <- video_pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(video_pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(video_pgls_WinFreq_Mass2_3)$coef[2,4]
text(25, 0.55, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(25, 0.55 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(25, 0.55 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
dev.off()
Fit Lanchester’s laws as well as the body size-only and group size-only models to the live-only and video-only datasets:
live_group <- sheet4[sheet4$Group.A != sheet4$Group.B,]
video_group <- scored4[scored4$Group.A != scored4$Group.B,]
# Exclude cornetfish:
live_group <- live_group[live_group$Species.A != "Cornetfish" &
live_group$Species.B != "Cornetfish",]
# Number of data points:
nrow(live_group)
## [1] 59
nrow(video_group)
## [1] 24
# Get body size info and use it to test the four models:
lanchester_on_frames <- function(dataframe) {
n <- ncol(dataframe)
dataframe$Group.A <- as.numeric(dataframe$Group.A)
dataframe$Group.B <- as.numeric(dataframe$Group.B)
for(i in 1:nrow(dataframe)) {
dataframe[i, n+1] <- lengthmass$Mass[lengthmass$Species == as.character(dataframe$Species.A[i])]
dataframe[i, n+2] <- lengthmass$Mass[lengthmass$Species == as.character(dataframe$Species.B[i])]
dataframe[i, n+3] <- lengthmass$Mass2_3[lengthmass$Species == as.character(dataframe$Species.A[i])]
dataframe[i, n+4] <- lengthmass$Mass2_3[lengthmass$Species == as.character(dataframe$Species.B[i])]
if (dataframe$Group.A[i]*dataframe[i, n+1] > dataframe$Group.B[i]*dataframe[i, n+2]) {
dataframe[i, n+5] <- dataframe$Species.A[i]
} else {
dataframe[i, n+5] <- dataframe$Species.B[i]
}
if ((dataframe$Group.A[i]^2)*dataframe[i, n+1] > (dataframe$Group.B[i]^2)*dataframe[i, n+2]) {
dataframe[i, n+6] <- dataframe$Species.A[i]
} else {
dataframe[i, n+6] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+5])) == TRUE) {
dataframe[i, n+7] <- 1
} else {
dataframe[i, n+7] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+6])) == TRUE) {
dataframe[i, n+8] <- 1
} else {
dataframe[i, n+8] <- 0
}
if (dataframe$Group.A[i]*dataframe[i, n+3] > dataframe$Group.B[i]*dataframe[i, n+4]) {
dataframe[i, n+9] <- dataframe$Species.A[i]
} else {
dataframe[i, n+9] <- dataframe$Species.B[i]
}
if ((dataframe$Group.A[i]^2)*dataframe[i, n+3] > (dataframe$Group.B[i]^2)*dataframe[i, n+4]) {
dataframe[i, n+10] <- dataframe$Species.A[i]
} else {
dataframe[i, n+10] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+9])) == TRUE) {
dataframe[i, n+11] <- 1
} else {
dataframe[i, n+11] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+10])) == TRUE) {
dataframe[i, n+12] <- 1
} else {
dataframe[i, n+12] <- 0
}
if (dataframe$Group.A[i] > dataframe$Group.B[i]) {
dataframe[i, n+13] <- dataframe$Species.A[i]
} else {
dataframe[i, n+13] <- dataframe$Species.B[i]
}
if (dataframe[i, n+1] > dataframe[i, n+2]) {
dataframe[i, n+14] <- dataframe$Species.A[i]
} else {
dataframe[i, n+14] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+13])) == TRUE) {
dataframe[i, n+15] <- 1
} else {
dataframe[i, n+15] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+14])) == TRUE) {
dataframe[i, n+16] <- 1
} else {
dataframe[i, n+16] <- 0
}
}
colnames(dataframe)[(n+1):(n+16)] <- c("BodySize.A", "BodySize.B", "BodySize2_3.A",
"BodySize2_3.B", "Linear.Law.Prediction",
"Square.Law.Prediction", "Linear.True",
"Square.True", "Linear.Law2_3.Prediction",
"Square.Law2_3.Prediction", "Linear2_3.True",
"Square2_3.True", "Group.Size.Only.Prediction",
"Body.Size.Only.Prediction",
"Group.Size.Only.True", "Body.Size.Only.True")
return(dataframe)
}
live_group <- lanchester_on_frames(live_group)
video_group <- lanchester_on_frames(video_group)
# Proportion of successful predictions:
sum(live_group$Linear.True)/nrow(live_group)
## [1] 0.6949153
sum(live_group$Square.True)/nrow(live_group)
## [1] 0.5254237
sum(live_group$Linear2_3.True)/nrow(live_group)
## [1] 0.6101695
sum(live_group$Square2_3.True)/nrow(live_group)
## [1] 0.4237288
sum(live_group$Group.Size.Only.True)/nrow(live_group)
## [1] 0.2542373
sum(live_group$Body.Size.Only.True)/nrow(live_group)
## [1] 0.6779661
sum(video_group$Linear.True)/nrow(video_group)
## [1] 0.375
sum(video_group$Square.True)/nrow(video_group)
## [1] 0.08333333
sum(video_group$Linear2_3.True)/nrow(video_group)
## [1] 0.1666667
sum(video_group$Square2_3.True)/nrow(video_group)
## [1] 0.04166667
sum(video_group$Group.Size.Only.True)/nrow(video_group)
## [1] 0
sum(video_group$Body.Size.Only.True)/nrow(video_group)
## [1] 0.75
# Test significance:
dbinom(sum(live_group$Linear.True), size = nrow(live_group), prob = 1/2)
## [1] 0.001123269
dbinom(sum(live_group$Square.True), size = nrow(live_group), prob = 1/2)
## [1] 0.09596023
dbinom(sum(live_group$Linear2_3.True), size = nrow(live_group), prob = 1/2)
## [1] 0.02501637
dbinom(sum(live_group$Square2_3.True), size = nrow(live_group), prob = 1/2)
## [1] 0.05253438
dbinom(sum(live_group$Group.Size.Only.True), size = nrow(live_group), prob = 1/2)
## [1] 6.920778e-05
dbinom(sum(live_group$Body.Size.Only.True), size = nrow(live_group), prob = 1/2)
## [1] 0.002423897
dbinom(sum(video_group$Linear.True), size = nrow(video_group), prob = 1/2)
## [1] 0.07793331
dbinom(sum(video_group$Square.True), size = nrow(video_group), prob = 1/2)
## [1] 1.645088e-05
dbinom(sum(video_group$Linear2_3.True), size = nrow(video_group), prob = 1/2)
## [1] 0.000633359
dbinom(sum(video_group$Square2_3.True), size = nrow(video_group), prob = 1/2)
## [1] 1.430511e-06
dbinom(sum(video_group$Group.Size.Only.True), size = nrow(video_group), prob = 1/2)
## [1] 5.960464e-08
dbinom(sum(video_group$Body.Size.Only.True), size = nrow(video_group), prob = 1/2)
## [1] 0.008022547
How many interactions involving damselfish do we have in the live-only dataset and how many in the video-only dataset?
nrow(live_group[live_group$Species.A == "Damselfish" | live_group$Species.B == "Damselfish",])/nrow(live_group)
## [1] 0.3559322
nrow(video_group[video_group$Species.A == "Damselfish" | video_group$Species.B == "Damselfish",])/nrow(video_group)
## [1] 0.7916667